TECHNICAL  REPORT 


Defense  Threat  Reduction  Agency 
8725  John  J.  Kingman  Road,  MS-6201 
Fort  Belvoir,  VA  22060-6201 


DTRA-TR-14-49 


Meso-Scale  Experimental  &  Numerical 
Studies  for  Predicting  Macro-scale 
Performance  of  Advanced  Reactive 
Materials  (ARMs) 

Distribution  Statement  A.  Approved  for  public  release;  distribution  is 
unlimited. 


April  2015 


HDTRA1 -10-1 -0038 


Naresh  Thadhani  et  al. 

Prepared  by: 

Georgia  Tech  Research  Corporation 

School  of  Materials  Science  and 

Engineering 

771  Ferst  Drive 

Atlanta,  GA  30332 


DESTRUCTION  NOTICE 


Destroy  this  report  when  it  is  no  longer  needed. 
Do  not  return  to  sender. 


PLEASE  NOTIFY  THE  DEFENSE  THREAT  REDUCTION 
AGENCY,  ATTN:  DTRIAC/  J9STT,  8725  JOHN  J.  KINGMAN  ROAD, 
MS-6201,  FT  BELVOIR,  VA  22060-6201,  IF  YOUR  ADDRESS 
IS  INCORRECT,  IF  YOU  WISH  THAT  IT  BE  DELETED  FROM  THE 
DISTRIBUTION  LIST,  OR  IF  THE  ADDRESSEE  IS  NO 
LONGER  EMPLOYED  BY  YOUR  ORGANIZATION. 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Suhithi  Peiris 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

SAR 

54 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

703-767-4732 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


MESO-SCALE  EXPERIMENTAL  &  NUMERICAL  STUDIES  FOR  PREDICTING  MACRO¬ 
SCALE  PERFORMANCE  OF  ADVANCED  REACTIVE  MATERIALS  (ARMS) 


Grant/Award#:  HDTRA1-10-1-0038 

Jennifer  Breidenich,  Arun  Gokhale,  Manny  Gonzales,  Ashok  Gurumurthy,  and 

Naresh  Thadhani 

School  of  Materials  Science  and  Engineering,  Georgia  Institute  of  Technology 

and 


Jason  Quenneville 

Spectral  Sciences  Inc.,  Spectral  Sciences  Inc.,  Burlington,MA  01803 


P.l.  Name:  Naresh  Thadhani 
Organization:  Georgia  Tech  Research  Corporation 
Address:  School  of  Materials  Science  and  Engineering,  Georgia  Institute  of 

Technology,  771  Ferst  Drive,  Atlanta,  GA  30332-0245 


Final  Report 
April  28,  2014 


2 


TABLE  OF  CONTENTS 


ABSTRACT  4 

1.  INTRODUCTION  4 

2.  BACKGROUND  5 

2.1  Characteristics  of  Anaerobic  Reactions  5 

2.2  Configurational  Effects  on  Anaerobic  Chemical  Reactions  5 

2.3  Characteristic  Scales  for  Defining  Meso-Scale  Processes  6 

2.4  Microstructure  Generation  and  Quantification  7 

3.  TECHNICAL  APPROACH  AND  SIGNIFICANT  RESULTS  9 

3.1  Impact  Initiated  Combustion  of  Aluminum  10 

3.2  Impact-Initiated  Reactions  in  Ti+AI+B  Mixtures  Under  Uniaxial  Stress  16 

3.3  Shock-Compression  of  Ti+2B  Mixtures:  Experiments  and  Simulations  19 

4.  Statistical  Representation,  Reconstruction,  Simulations  of 

Meso-Scale  3D  Microstructures  of  Reactive  Materials  25 

4.1  Characteristics  of  Real  Microstructures  29 

4.2  The  Simulation  Method  30 

4.3  Parameterization  of  Shape  32 

4.4  Parameterization  of  Clustering  33 

4.5  DoE  for  Delineating  Effects  of  Microstructural  Parameters  34 

4.6  Numerical  Proof  of  Anisotropy  of  Compact  Microstructures  35 

5.  Non-equilibrium  MD  Simulations  of  Shocked  Ni/Al  Nanolaminates  36 

5.1  Simulation  of  Ni+Al  Nanolaminates  36 

5.2  Toward  MD  Simulation  of  Shock-Induced  Reactions  of  Ti  and  B  43 

BIBLIOGRAPHY  50 


3 


Grant/ Award  ft: 
PI  Name: 
Organization: 
Project  Title: 


HDTRA1-10-1-0038 

Naresh  Thadhani 

Georgia  Institute  of  Technology 

Meso-scale  Experimental  &  Numerical  Studies  for  Predicting  Macro-scale 
Performance  of  Advanced  Reactive  Materials  (ARMs) 


ABSTRACT 

Our  collaborative,  multi-disciplinary  and  multi-investigator  research  project,  combines 
numerical  simulations  with  time-resolved  impact  experiments,  to  determine  the  meso-scale 
mechanisms  of  impact-initiated  reactions  in  Ti-AI-B  based  reactive  materials  in  the  form  of 
compacts  of  powders  of  different  sizes  and  morphologies.  The  major  goal  is  to  delineate  how 
processes  of  localized  deformation  and  flow,  or  fracture  and  fragmentation,  create  the 
dispersion  and  mixing  of  reactants  and  provide  the  thermal  and  mass  transport  resulting  in 
highly  exothermic  Ti  and  B  intermetallic  reaction  and/or  oxidation  reaction  (Al  or  B).  The  time 
sequence  of  events  associated  with  reaction  initiation  as  well  as  subsequent  interactions  with 
the  external  environment  are  also  being  determined.  More  specifically,  the  influence  of 
material-inherent  elastic/plastic  properties  and  reactant  configuration  (e.g.,  porosity, 
morphology,  spacing,  distribution,  etc.,  created  from  a  library  of  microstructures)  on  meso- 
scale  mechanistic  processes  and  the  resulting  macro-scale  performance  is  being  established. 
The  combination  of  experiments  and  numerical  simulations  that  capture  morphologies  of  real 
materials  through  the  creation  of  synthetic  microstructures  from  a  library  of  constituent 
morphologies  is  used  to  correlate  processes  across  various  length  scales  and  predict 
fundamental  spatial  and  temporal  understanding  of  meso-scale  intermetallic  and/or  oxidation 
reaction  processes.  The  fundamental  predictive  studies  have  the  potential  to  enable  designs  of 
the  meso-scale  structure  of  advanced  reactive  materials. 


1.  INTRODUCTION 

The  research  performed  in  this  work  was  aimed  at  generating  the  fundamental  spatial 
and  temporal  understanding  of  the  meso-scale  mechanisms  of  impact-initiated  anaerobic 
chemical  reactions  in  intermetallic-based  reactive  materials  in  the  form  of  nano/micro-sized 
particulate  compacts  and  layered  structures.  The  specific  objectives  of  the  proposed  work  are 
to  (a)  delineate  how  processes  of  localized  deformation  and  flow,  or  fracture  and 
fragmentation,  create  the  dispersion  and  mixing  of  reactants  and  provide  the  thermal  and  mass 
transport  needed  for  anaerobic  chemical  reaction  initiation;  (b)  determine  the  time  sequence 
of  events  associated  with  reaction  initiation,  as  well  as  subsequent  interactions  with  external 
environment;  and  (c)  establish  the  role  of  material-inherent  elastic/plastic  properties  and 
reactant  configuration  (e.g.,  porosity,  morphology,  distribution,  etc.)  on  the  mechanisms  and 
kinetics  of  meso-scale  processes  and  the  resulting  macro-scale  energetic  performance  of  such 
anaerobic  reactions.  These  objectives  have  been  pursued  through  four  distinct  activities, 
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integrated  to  generate  a  predictive  understanding  of  the  meso-scale  mechanisms  of  impact- 
initiated  reactions  in  Ti-Al-B  based  reactive  materials. 

2.  BACKGROUND 

2.1  Characteristics  of  Anaerobic  Reactions 

The  energy  release  characteristics  of  intermetallic-forming  reactive  materials  combined 
with  their  mechanical  properties  makes  them  highly  desirable  in  defense  weapons 
applications.1  Reactive  materials  are  inherently  heterogeneous,  with  meso-scale  variations  in 
nano-  to  micro-meter  size  of  reactants  of  various  forms/morphologies.  Shock  initiation  of 
reactions  is  dominated  by  physical,  mechanical,  and  chemical  processes  that  occur  in  the 
nanosecond  time-scale  of  the  rise  to  the  peak  pressure,  to  the  tens  of  microsecond  time  scale 
of  temperature  equilibration.2  Understanding  the  influence  of  meso-scale  spatial  variations  and 
probing  processes  in  real  time  is  essential  for  predicting  the  macro-scale  reaction  response. 
Time-resolved  experiments  employing  conventional  diagnostics  provide  nanosecond  temporal 
resolution,  yet  there  is  a  limit  to  the  type  of  information  obtained.  Hence,  observations  of 
indirect  inference  of  reaction  have  often  been  utilized  based  on  expanded  P-V  compressibility 
states3  or  increased  shock  speeds.4 

Details  of  meso-scale  processes  in  reactant  particles  are  often  investigated  through 
computations.  Discrete-element  computations5  have  provided  an  excellent  parametric  study  of 
the  influence  of  internal  void  content,  strain-rate,  relative  yield  strength,  and  melting  on  the 
degree  of  mixing  and  reactions  in  regions  of  highly  localized  deformation  (or  shear-bands).  The 
results  have  revealed  an  optimum  strength  ratio  between  components  of  most  effective  inter¬ 
reactant  mixing,  the  extent  of  which  increases  with  shear  strain-rate.  Melting  occurring  more 
rapidly  at  high  strain-rates,  leads  to  loss  of  strength,  and  increases  strength  ratio,  producing  a 
negative  effect  on  mass  mixing  such  that  the  chemical  reaction  is  hindered.  Computations6  with 
simulated  (2-D)  microstructures  have  been  used  to  describe  the  evolution  of  reactant  and 
product  configurations  following  shock  compression.  Computations  with  3-D  simulated 
microstructures  have  also  been  performed.7  However,  similar  to  the  discrete  element  method, 
these  computations  have  employed  idealized  reactant  geometries  (spherical,  cubic,  angular) 
and  their  distributions,  which  are  disadvantaged  in  the  case  where  the  characteristics  of 
heterogeneities  are  far  from  ideal,  and  where  variations  exist  across  different  length  scales, 
particularly  if  the  goal  is  to  perform  predictive  modeling  of  macro-scale  performance. 

2.2  Configurational  Effects  on  Anaerobic  Chemical  Reactions 

Intermetallic-based  reactive  materials  are  non-explosive  and  non-detonable,  but  can 
rapidly  release  large  amounts  of  the  anaerobic  heat  of  reaction.8  While  processes  of  thermal 
ignition  in  such  reactive  materials  are  well  understood,9  the  mechanisms  of  impact-initiated 
reactions  and  their  kinetics  remain  to  be  established.  Our  prior  work10,11,12  showed  that  shock- 
initiation  of  chemical  reactions  is  influenced  by  differences  in  the  elastic/plastic  properties  of 
reactant  materials  and  their  morphological  characteristics.  Under  uniaxial-strain  loading,  micro- 
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size  spherical  Ni+Al  powder  mixtures  showed  crush-up  to  full  density  at  ~0.5  GPa,  while  the 
crush-up  stress  for  flake-Ni+Al  mixtures  of  similar  size  is  ~2  GPa,  and  that  of  nano-sized  Ni+Al 
was  ~6  GPa.13  However,  flake-Ni+Al  mixtures  reacted  at  the  lowest  threshold  stress,  due  to  the 
additional  bending  and  buckling  modes  that  cause  the  flake-Ni  particles  to  undergo  the  same 
degree  of  plastic  deformation  as  the  lower  strength  Al  particles,  thereby  promoting  mixing  and 
consequent  reaction  initiation.  Likewise,  Ti  and  Si  powder  mixtures  have  shown  variable 
reactivity  as  a  function  of  size,  dictated  by  the  ability  of  the  covalently  bonded  Si  particles  to 
undergo  deformation  or  fracture.8  Barring  the  fine  (<2  pm)  powders  which  tend  to 
agglomerate,  Si  particles  of  10-20  pm  size  undergo  plastic  deformation  simultaneously  with  the 
Ti  particles  of  the  same  size,  while  coarse  Si  particles  (>100  pm)  fracture  and  become  entrapped 
within  the  deformed  Ti.  Lack  of  flow  of  Si  particles  due  to  agglomeration,  or  fracture  and 
entrapment,  restricts  mixing  and  subsequent  reaction,  while  simultaneous  deformation  and 
flow  of  Si  with  Ti,  promotes  mixing  and  subsequent  mechanochemical  reaction  at  a  threshold 
shock-pressure  of  ~1  GPa.  A  threshold  shock-energyu  criterion  for  reaction  initiation  has  also 
been  postulated  for  Ti+Si  mixtures  with  varying  porosity  levels.  The  energy  threshold 
corresponds  to  the  work  of  shock  compression  leading  to  a  mean  equilibrated  (bulk) 
temperature  rise  exceeding  that  needed  for  thermally-initiated  reaction.  However,  the 
temperature  equilibration  time-scale  (tens  of  microseconds  to  milliseconds)  significantly 
exceeds  the  shock-pulse  duration.  Hence,  the  threshold  shock-energy  criterion  is  applicable  to 
shock-assisted  thermochemical  reactions  occurring  after  unloading,  and  contrasting  with  shock- 
induced  mechanochemical  reactions  that  occur  during  the  shock  state.2 

The  varied  mechanical  and  chemical  responses  described  above  are  controlled  by  meso- 
scale  processes,  which  in  turn  are  influenced  by  (a)  stress  states  and  pressure/temperature 
histories  generated  due  to  the  macro-scale  geometry,  impact  configuration,  and  loading 
conditions,  (b)  variations  in  inter-reactant  contact  areas  and  local  stress/strain  distributions 
caused  by  morphological  form  (shape,  size,  distribution)  of  reactants,  and  (c)  inherent 
differences  in  elastic  and  plastic  properties  of  reactants.  Probing  of  meso-scale  processes  is 
therefore,  essential  for  developing  a  complete  understanding  of  the  mechanisms  and  kinetics  of 
anaerobic  reactions  and  the  energy  release  characteristics  of  reactive  materials.  The 
understanding  generated  can  allow  correlation  of  the  macro-scale  performance  (in  terms  of  the 
extent  of  reaction  and  associated  temperature  and  pressure  changes)  with  the  structure,  and 
also  make  it  possible  to  reverse  design  the  reactants  architecture  for  desired  macro-scale 
properties  and  performance  characteristics. 

2.3  Characteristic  Scales  for  Defining  Meso-scale  Processes 

Probing  meso-scale  processes  in  a  macroscopic  bulk  reactive  material  requires  defining  a 
representative  volume  element  (RVE)  that  it  is  large  enough  to  capture  the  spatial  scales, 
distributions,  and  morphological  features  of  reactants.  The  RVE  however,  can  be  quite  different 
for  nanometer  or  micrometer  size  reactant  particles/grains/layers/laminates.  As  illustrated  in 
Fig.  1,  the  description  of  meso-scale  processes  for  an  RVE  also  requires  defining  the  stress-state 
and  strain-rate  bounds  from  the  higher  level  associated  with  macro-scale  sample 
geometry/configuration,  and  likewise  the  lower  level  local  stress-strain  distributions,  processes 
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of  plastic  flow  and/or  fracture,  mixing/transport  of  reactants 
occurring  at  the  particle-scale,  as  well  as  the  interface-scale. 
Use  of  real  imported  microstructures  enables  such  a  linkage  by 
providing  a  realistic  representation  of  the  size  and  form  of  the 
constituents.  It  also  makes  it  possible  to  probe  and  correlate 
processes  across  the  length  scales.  The  use  of  imported 
microstructures  to  describe  the  response  of  materials  has  been 
used  to  address  a  wide  variety  of  materials  applications, 
including  predictions  of:  (a)  3-D  fluid  flow  and  permeability  in 
asphalt  pavements,15  (b)  thermal  conductivity  of  thermal  barrier 
coatings,16  (c)  dynamics  of  fracture  path  and  crack  propagation 
in  Al+SiC  composites17  and  TiB2+Al203  two-phase  ceramics;18  (d) 
mechanics  of  biological  materials  for  simulating  brain  tumor 
growth  and  its  impact  on  surrounding  tissues;19  and  (e)  linear  elasticity  and  thermal 
conductivity  of  two-phase  polycrystalline  ice  cream20  to  design  the  structure  needed  for 
prescribed  taste  and  texture.  Shock-compression  response  of  materials  has  been  studied  by 
Benson,21  who  utilized  a  pixel-by-pixel  correlation  technique  to  import  the  microstructure  of 
HMX  into  Raven,  22  to  simulate  local  temperature  variations  during  shock-propagation  in 
faceted  and  equiaxed  HMX  crystals.  Benson23  has  also  used  real  microstructures  to  accurately 
reproduce  the  equation  of  state  of  Cu.  Use  of  real  microstructures  has  also  been  made  to 
describe  the  shock-compression  response  of  reactive  mixtures  under  uniaxial-strain24,25  and 
uniaxial-stress26  loading.  These  studies  have  demonstrated  that  real  imported  microstructures 
allow  more  accurate  predictions  in  systems  with  inherent  heterogeneities  that  significantly 
influence  the  phenomena  or  processes  being  investigated.  Use  of  real  microstructures  also 
makes  it  possible  to  correlate  mathematical  representations  of  the  microstructure  for 
correlation  with  performance,  thereby  attaining  the  most  robust  predictive  capability  and 
enabling  the  reverse  design  of  meso-scale  structure  for  required  performance  characteristics, 
which  will  be  the  focus  of  the  proposed  work. 


2.4  Microstructure  Generation  and  Quantification 

Meso-scale  microstructures  of  multi-phase  Ti+AI+B  intermetallic-forming  reactive 
materials  are  three-dimensional  (3-D)  and  stochastic,  with  geometrical  shapes  that  are  complex, 
morphological  orientations  that  are  not  necessarily  isotropic,  distributions  of  reactants  that  are 
not  necessarily  uniform-random,  spatial  correlations  and  clustering  that  exist  frequently,  and 
interactions/correlations  among  features  occurring  at  different  length  scales.  The  geometric 
complexities  present  challenges  for  quantitative  characterization,  statistical  representation,  and 
simulations  of  microstructures,  but  also  provide  opportunities  to  optimize  materials  properties 
and  enhance  their  performance.  In  this  work,  we  have  incorporated  real  reconstructed  high- 
resolution  large-volume  3D  microstructures  as  RVEs  to  capture  the  complex  meso-scale 
geometries  of  reactants.  It  build  on  our  prior  work  on  montage  serial  sectioning27,28,29  and 

qn  3^  33  3  j  qr 

reconstructing  images  '  '  for  FE  simulations,  '  '  to  establish  a  framework  for  predicting 
the  macro-scale  performance  of  advanced  reactive  materials. 
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Statistical  Representation  of  Microstructures:  A  microstructure  is  an  ensemble  of 
volumes,  internal-surfaces,  lines,  and  points.  Depending  on  its  dimensionality,  each 
microstructural  feature  (particle,  pore,  grain)  has  an  associated  size,  shape,  volume,  surface 
area,  curvature,  orientation,  and  location.  Specification  and  statistical  representation  of  the 
ensemble  geometry  requires  quantification  of  volume  fraction,  total  surface  area,  size 
distributions,  etc.,  as  well  as  descriptors  of  spatial  arrangements,  morphological  anisotropy,  and 
topology.  Stereology  techniques  are  well  established  for  efficient/unbiased  estimation  of 
volume  fraction,  total  surface  area  per  unit  volume,  length  per  unit  volume,  number  per  unit 
volume,  and  integral  mean  curvature.36,37  Statistical  distribution  functions  such  as  n-point 
correlation  functions  and  lineal  path  probability  distributions  go  a  step  further  and  provide 
information  on  spatial  arrangements  (clustering,  etc.)  and  anisotropy,  for  statistical 
representation  of  microstructure38,39  and  for  use  as  robust  descriptors  for  computer  simulations 
of  realistic  microstructures.  For  a  3-D  microstructure  an  n-point  correlation  function  PKjt  ^  Jrv 
r2,..r„)  is  the  probability  that  a  randomly  located  polyhedron  having  n-vertices,  edge  lengths  rj, 
r 2/  ...,  and  specific  orientation  superimposed  on  a  microstructure  is  such  that  its  first  vertex  lies 
in  the  phase  /',  second  vertex  lies  in  the  phase  j,  and  so  on.40  Figure  2  illustrates  examples 
showing  (a)  two-point,  (b)  three-point,  and  (c)  four-point  correlation  functions  in  Al+SiC  two- 
phase  material.  Another  statistical  descriptor  is  the  lineal  path  probability  distribution  /.«  (r,  9, 
(p),  which  is  the  probability  that  a  randomly  located  line  in  the  microstructural  domain  of 
interest  having  length  rand  angular  orientation  (9,  <p)  in  the  3-D  space  is  completely  contained 
in  the  phase  /.  The  variation  of  this  probability  with  line  length  r  and  its  angular  orientation  (9, 
(p)  contains  microstructural  information  not  available  with  correlation  functions.  In  the 
proposed  research,  unbiased  stereology  and  image  analysis  based  computation  techniques42,43 
will  be  used  to  obtain  first-order  global  microstructural  properties  (volume  fraction,  surface 
area  per  unit  volume, 
etc.),  correlation 

functions,  and  lineal 
path  probability 

distributions  for 

detailed  statistical 
representation  of  the 
microstructures,  and  to 
establish  quantitative 
relations  between 
processing,  structure,  and  properties. 

Reconstruction  of  3-D  Microstructures:  Random  independent  2-D  metallographic 
sections  do  not  contain  all  information  concerning  the  true  3-D  microstructural  geometry. 
Information  regarding  true  3-D  particle/feature  shapes/morphologies,  spatial  clustering  and 
correlations,  and  topological  aspects  such  as  3-D  connectivity  (e.g.,  of  pores)  does  not  reside  in 
random  independent  2-D  metallographic  sections.40  Hence,  incorporation  of  high-resolution 
large  segments  of  realistic  3-D  microstructural  images  as  RVEs  in  computations,  is  required  to 
account  for  effects  of  these  important  attributes  on  material  behavior.33  We  have  developed  an 
efficient  montage-based  serial  sectioning  technique41  for  generation  of  large  volume  segments 


SI/ 


Fig.  2  Images  showing  (a)  two-,  (b)  three-,  (c)  four-point  correlation  functions 
corresponding  to  vertices  in  two-phase  material  SiC  particles  in  A1  matrix.3* 
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of  high  resolution  3-D  microstructures  from  stacks  of  aligned  montage  serial  sections,  and 
successfully  applied  to  reconstruct  3-D  microstructures  of  numerous  materials,25’33'34  which  will 
be  used  in  the  proposed  work.  Figure  3  illustrates  an  example  of  3-D  microstructures 
reconstructed  from  montage  serial  sections.  Complete  reconstructed  3-D  volumes  are  much 
larger  than  those  depicted.  High-resolution  large  volume  segments  of  3-D  meso-scale 
microstructures  of  reactive  materials  will  be  constructed  using  this  reconstruction  methodology 


Fig.  3.  Reconstruction  of  3D  microstructure  of  Al-alloy  composite  containing  SiC  particles  having  8. 1  particle 
size  ratio  (PSR).  (a)  Stack  of  small  number  of  digitally  compressed  aligned  montage  serial  sections  created  by 
stitching  100  contiguous  microstructural  fields,  (b).  Each  field  of  view  hi  (a)  was  acquired  at  resolution  depicted 
in  (b);  (c)  Example  of  small  segment  of  reconstructed  3-D  microstructure  depicting  clustering  of  SiC  particles.41 


3.  TECHNICAL  APPROACH  AND  SIGNIFICANT  RESULTS 

Advanced  reactive  materials  based  on  the  intermetallic-forming  Ti+B+AI  system  were 
investigated  in  this  work.  The  mixtures  of  the  reactive  elemental  systems  represent 
constituents  that  undergo  highly  exothermic  anaerobic  reactions,  while  also  making  it  possible 
for  residual  reactants  to  become  heated  for  subsequent  (combustion)  reactions.  Various 
possible  intermetallic  and  elemental  combustion  reactions  can  occur  during  impact  loading  of 
Ti+B+AI  mixtures.  The  equilibrium  heat  of  reaction10  of  combustion  of  Al  (AHR=-7422  cal/g), 
Ti(AHR=-4714  cal/g),  and  B(AHR=-14050  cal/g)  is  higher  than  that  of  intermetallic  reactions,  e.g., 
AI+2B  (AHr=-371  cal/g),  2AI+Ti(AHR=-314  cal/g),  and  2B+Ti(AHR=-1320  cal/g).  The  threshold  for 
reaction  initiation  is  more  dependent  on  the  volumetric  distribution  of  reactants,  their  intrinsic 
properties,  and  morphologies,  rather  than  the  heat  of  reaction.  The  constituents  have  varying 
intrinsic  elastic/plastic  properties,  and  the  differences  in  these  properties  affect  the  extent  of 
their  deformation  and  fracture,  the  degree  of  their  dispersion  and  mixing,  and  resulting 
threshold  for  reaction  initiation.  In  the  case  of  Ti+AI+B,  the  large  difference  in  flow  stress 
between  Ti  and  Al  (Aay>100  MPa)  can  lead  to  extensive  flow  of  Al  without  significant  straining 
of  Ti,  which  limits  their  mixing.  However,  the  brittle  behavior  and  high  sound  speed  of  B 
(ACo~10km/s),  in  spite  of  its  low  strength,  can  force  its  dispersion  and  mixing  with  Ti.42 

In  addition  to  the  constituent  characteristics,  the  onset  of  reaction  is  also  influenced  by 
the  state  of  loading  involving  uniaxial  stress  or  uniaxial  strain  conditions.  Hence,  the  overall 
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approach  adopted  in  this  work  was  to  employ  controlled  gas  gun  experiments  investigate  (a) 
impact-initiated  combustion  reaction  of  aluminum  powders  of  varying  morphology  under 
conditions  of  uniaxial-stress  loading,  (b)  impact-initiated  reactions  in  Ti+Al+B  mixtures  (Ti+2B 
mixed  with  varying  contents  of  Al)  under  conditions  of  uniaxial  stress  loading,  and  (c)  shock- 
compression  of  binary  Ti+2B  powder  mixtures  under  conditions  of  uniaxial-strain  loading. 
Meso-scale  simulations  with  actual  (real)  and  simulated  microstructures  were  also  performed 
with  synthetic  (simulated)  microstructures  generated  from  a  library  of  constituent 
morphologies  of  particles  and  validation  through  two-point  nearest  neighbor  correlations.  The 
significant  results  obtained  during  the  course  of  this  work  are  detailed  next. 


3.1  Impact-initiated  combustion  of  aluminum 

The  mechanisms  of  impact  initiated  combustion  of  aluminum  powders  of  different  sizes, 
morphologies,  and  pre-conditioning  were  investigated  using  uniaxial-stress  rod-on-anvil  impact 
loading  with  the  7.62-mm  single-stage  gas-gun.  Aluminum  powders  of  a  range  of  particle  sizes 
and  exposed  to  mechanical  pre-activation  were  studied  to  analyze  the  aluminum  combustion 
reaction  occurring  during  crush-up  or  subsequent  compact  deformation.  The  initiation 
mechanism,  either  thermochemical,  due  to  densification  and  deformation  of  the  compact,  or 
mechanochemical,  due  to  particle-particle  interactions,  or  a  combination  of  both  mechanisms, 
were  examined.  A  clear  understanding  of  the  mechanisms  responsible  for  the  combustion 
reaction  in  aluminum  powder  compacts  can  aid  in  the  determination  of  the  primary  reaction 
process  during  the  reactions  of  binary  and  ternary  Ti+Al+B  intermetallic-forming  mixtures. 

Reactivity  of  Aluminum  Powders  of  Different  Sizes  and  Size  Distributions:  The  initiation  of  the 
aluminum  combustion  reaction  under  high-strain-rate  impact-loading  conditions  was 
investigated  by  determining  the  onset  of  reaction  as  a  function  of  aluminum  particle  size, 
particle  size  distribution,  density  of  the  compacts,  and  impact  velocity.  Spherical,  as-received 
aluminum  powders  of  four  different  diameters  were  investigated,  a  ~3pm  diameter  powder 
(Valimet  H-2,  Figure  4a),  a  ~9pm  diameter  powder  (Valimet  H-5,  Figure  4b),  a  ~30p.m  aluminum 
powder  (Valimet  H-15,  Figure  4c),  and  a  ~75pm  aluminum  powder  (Valimet  H-60,  Figure  4d).  A 
Saltykov  analysis  of  the  particle  size  distribution  gave  a  normal  particle  size  distribution  for  the 
as-received  powder  sets.  Additionally,  mixtures  with  a  bimodal  size  distribution  were  created 
by  combining  H-60  and  H-2  powders  in  a  1:3  ratio  (Figure  4e)  and  a  3:1  ratio  (Figure  4f)  and 
mixing  in  a  V-blender. 
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Figure  4.  Images  of  aluminum  powder  mixtures  tested.  SEM  images  of  Valimet  H-2  (a),  H-5  (b),  H-15  (c),  and  H-60 
(d).  Light  microscope  images  of  bimodal  mixtures  of  H-60  and  H-2, 1:3  H-60:H-2  (e)  and  3:1  H-60:H-2  (f). 


SEM  and  XRD  analysis  were  performed  on  both  the  as-received  and  highly-strained 
aluminum  powders.  XRD  line-broadening  analysis  employing  the  Williamson-Hall  method  was 
used  to  determine  the  retained  lattice  strain  in  each  type  of  impacted  sample. 

All  starting  powders  were  quasi-statically  pressed  into  compacts  of  2.0  ±  0.1mm  height 
and  3.19mm  diameter,  ranging  in  density  from  80%  to  99%  TMD.  The  surface  area  to  volume 
ratio  (Sv)  of  each  powder  type  was  measured  by  the  intersection  counting  method  from 
compacts  of  the  pressed  powders.  The  compacts  were  cross-sectioned,  polished,  etched,  and 
imaged  by  light  microscopy  (Leica  DM  IRM).  Cycloids  were  super-imposed  on  the  image,  and 
the  average  number  of  intersections  between  the  cycloid  and  the  microstructural  boundaries 
per  unit  of  the  cycloid  length  (PL)  were  related  to  the  surface  area  per  unit  volume,  Sv,  via  the 
relationship,  SV=2PL. 

The  aluminum  powder  compacts  were  affixed  to  the  lapped  end  of  a  copper  projectile 
(38.1mm  high  x  7.62mm  diameter)  via  24  hour  Hysol  and  impacted  against  a  hardened  steel 
anvil  using  the  7.62mm  gas  gun.  The  impact  event  was  imaged  via  a  high-speed  framing  camera 
(IMACON  200)  in  order  to  determine  the  initiation  of  combustion  reaction  during  impact. 

The  observation  of  light  in  several  consecutive  frames  captured  by  the  high-speed 
framing  camera  was  used  as  evidence  of  the  occurrence  of  impact-induced  combustion 
reaction,  as  shown  by  the  example  of  images  of  a  HEBM  aluminum  powder  compact,  illustrated 
in  Figure  5.  Each  powder  compact  was  tested  at  a  range  of  velocities,  from  140  m/s  to  320  m/s, 
to  determine  the  minimum  velocity  necessary  for  reaction  to  initiate  under  each  condition.  All 
experiments  were  performed  under  ambient  environmental  conditions. 
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Figure  5.  High-speed  images  of  a  HEBM 
aluminum  powder  compact  affixed  to  a  copper 
projectile  impacting  a  hardened  steel  anvil  at  302 
m/s.  Top  image  taken  immediately  prior  to 
impact.  Images  after  3,  6,  and  8  ns  following 
impact  show  light  emission  during  compaction 
and  deformation  of  pellet,  the  image  at  11  ns 
following  impact  show  light  emission  after 
complete  deformation  of  the  compact,  and  the 
final  image  at  1 5ns  after  impact  shows  sustained 
reaction  of  the  HEBM  pellet  as  the  Cu  projectile 
also  deforms  upon  making  contact  with  the  anvil. 


The  reactivity  of  each  powder  compact  type  was  plotted  against  the  measured  impact 
velocity  (Figure  6).  Reaction  occurrence  was  observed  to  vary  as  a  function  of  particle  size,  with 
lowest  velocities  necessary  for  reaction  initiation  ranging  from  190m/s  to  240m/s.  The  smallest 
particle  sizes,  H-2  (~3pm)  and  H-5  (~9pm),  showed  the  highest  initiation  threshold  velocity 
(~240m/s),  while  the  largest  diameter  particles,  H-60  (~75pm),  required  the  lowest  velocity 
(~190m/s)  to  initiate  reaction.  In  general,  reactivity  increased  with  particle  size.  The  1:3  and 
3:1  bimodal  mixtures  both  showed  reaction  at  velocities  (218m/s  and  223m/s,  respectively) 
between  the  largest  and  smallest  particle  sizes.  Additionally,  the  general  reactivity  trends  for 
each  powder  type  indicate  that  porosity  (or  compact  density)  also  plays  a  role  in  the  initiation 
of  the  combustion  reaction  in  the  Al  powder  compacts.  Compacts  with  densities  between  80 


and  90%  reacted  at  higher  velocities  than  those  with  densities  greater  than  90%  TMD. 
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Figure  6.  Impact  initiated  reactivity  of  aluminum  powder  compacts  as  a  function  of  particle  size,  particle  size 
distribution,  and  porosity  (compact  density). 
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Reactivity  of  Aluminum  exposed  to  mechanical  pre-activation:  The  effect  of  mechanical  pre¬ 
activation  on  the  impact-initiated  combustion  of  aluminum  powder  compacts  was  studied  by 
comparing  the  impact  response  of  as-received  aluminum  powders  to  those  strain-hardened  by 
either  ball-milling  or  by  high-strain  machining.  High  energy  ball  milled  (HEBM)  powders  were 
prepared  by  loading  5  g  of  Valimet  H-60  aluminum  powder  into  an  alumina  vial  and  ball  milling 
with  five  alumina  spheres  of  6.3mm  diameter,  using  a  SPEX  mill  for  10  minutes,  30  minutes,  and 
60  minutes  (Figure  7a-c).  During  the  milling  process,  the  spherical  powders  were  observed  to 
fracture,  plastically  deform,  and  cold  weld.  As  the  ball  milling  time  increased,  the  particles 
became  more 

flattened,  flake-like, 
and  gave  a  bimodal 
particle-size 
distribution  centered 
around  70  to  800pm 
diameter.  The  second 
mechanically  pre¬ 
activated  powder 
type,  highly  strained 
aluminum  6061  chips 
(figure  7d),  was 
prepared  by  Shankar 
et.  al.43  using  a  plane- 
strain  machining 
method.  These  chips 
contain  ultrafine 
grains  that  strengthen 
the  AI6061  and  are 
approximately  80pm  x  80pm  x  30pm  in  size.44 

The  impact-initiated  reactivity  of  each  powder  compact  type,  obtained  from  rod-on-anvil 
impact  experiments  described  above,  is  plotted  as  a  function  of  the  measured  impact  velocity, 
as  shown  in  Figure  8.  It  can  be  seen  that  after  processing  via  high  energy  ball  milling  for  10 
minutes,  the  reactivity  of  H-60  aluminum  decreases,  requiring  impact  velocities  of  at  least 
215m/s  to  react.  Yet,  as  ball  milling  times  are  increased  to  30  minutes,  the  reactivity  improves, 
with  reactions  observed  at  velocities  as  low  as  175m/s.  Compacts  of  Valimet  H-60  samples  ball- 
milled  for  60  minutes,  highly-strained  aluminum  platelets,  and  mixtures  of  H-60  and  platelet 
powders,  are  observed  to  be  less  reactive  than  H-60  powders,  with  reaction  velocity  thresholds 
in  the  range  of  207  to  213m/s. 
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Figure  8.  Impact  initiated  reactivity  of  aluminum  as  a  function  of  mechanical  pre-activation  and  porous  density. 


Possible  mechanisms  of  impact-initiated  combustion  of  aluminum  powder  compacts:  The 

minimum  (threshold)  impact  velocity  for  reaction  initiation  in  compacts  of  aluminum  powders 
is  found  to  correlate  with  particle  size,  mechanical  activation,  and  compact  density  based  on 
the  results  depicted  in  Table  I.  Mechanical  pre-activation  of  the  aluminum  powders  both  strain 
hardens  the  material  and  exposes  fresh  reactant  surfaces.  Strain  hardening  was  confirmed  from 
retained  elastic  strain  measured  using  the  Williamson-Hall  XRD  peak  broadening  method. 


Table  1.  Summary  of  aluminum 

powder  properties 

Particle  Type 

~Particle 

Diameter 

Range(pm) 

Surface  Area  per 
unit  Volume 

-i 

(mm  ) 

Lattice 

Strain 

Vickers 

Hardness 

Value 

Threshold 

Velocity 

(m/s) 

H-2 

1-7 

1885.7 

- 

- 

240 

H-5 

4-15 

824.3 

- 

- 

240 

H-15 

9-40 

441.9 

- 

- 

226 

H-60 

44-105 

105.7 

- 

190 

H-60  HEBM  lOmin 

44-250 

108.3 

.0013 

215 

H-60  HEBM  30min 

44-500 

166.9 

.0006 

175 

H-60  HEBM  60min 

60-750 

237.5 

.0005 

207 

Platelets 

80 

196.2 

.0049 

207 
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The  lattice  strain  for  highly-strained  aluminum  6061  platelets  was  significantly  greater 
than  the  lattice  strain  of  any  other  sample  tested,  but  it  resulted  in  decreased  reactivity,  in 
terms  of  a  higher  threshold  velocity.  On  the  other  hand,  Valimet  H-60  powders  exposed  to 
mechanical  pre-activation  via  HEBM,  showed  enhanced  reactivity  with  lower  threshold  velocity. 
Strain  hardening  by  mechanical  pre-activation  can  have  a  pronounced  effect  on  particle-particle 
interactions,  which  influences  the  energy  input  necessary  for  compaction  and  deformation  of 
the  powder  compact  during  impact.  Reactivity  of  the  highly  strained  aluminum  6061  powder 
compact  was  most  likely  decreased  because  the  harder,  stronger  particles  required  more 
energy  to  be  expended  during  compaction  and  deformation  of  the  higher  strength  particles, 
leaving  less  energy  available  for  reaction  initiation.  For  the  HEBM  H-60  particles,  high  energy 
ball  milling  for  a  short  period  of  time  causes  some  strain  hardening  as  observed  by  the 
measured  lattice  strain  but  also  breaks  the  particles  into  smaller  sizes  (up  to  30  minutes  milling) 
before  re-welding  and  agglomerating  the  powders  with  further  milling  time.  The  3-minute 
HEBM  samples  appear  to  exhibit  the  best  case  scenario  for  mechanical  pre-activation  resulting 
in  creation  of  fresh  reactant  surfaces  due  to  break-up  of  particles,  which  leaves  a  thinner,  more 
easily  penetrable  Al203  coating,  thereby  enhancing  the  reactivity  during  subsequent  impact. 

With  respect  to  particle  size,  the  threshold  velocity  for  reaction  of  compacts  of 
aluminum  powders  is  observed  to  decrease  with  increasing  particle  size.  The  1:3  bimodal 
particle  size  distribution  shows  decreased  reactivity  (higher  threshold  velocity).  Both  smaller 
particles  and  the  bimodal  particle  size  distribution  result  in  a  powder  compact  that  is  denser 
with  smaller  more  uniform  pore  spacing,  than  a  compact  composed  of  mostly  larger  particles. 
Compacts  containing  larger  particles  have  larger,  more  randomly  distributed  pores  and  are 
therefore  subjected  to  increased  plastic  deformation  and  heating  during  the  compaction  and 
deformation  process  during  impact.  Although  more  energy  is  necessary  to  plastically  deform 
the  particles,  the  enhanced  heating  serves  to  augment  the  reaction  in  the  larger  particles. 

The  validity  of  these  possible  mechanisms  for  impact  initiation  of  combustion  reactions 
in  the  various  aluminum  powder  compacts  can  be  further  verified  by  employing  a  transparent 
sapphire  anvil  (Figure  9),  instead  of  the  high  strength  steel  anvil.  Use  of  a  transparent  sapphire 
anvil  allows  observation  of  reaction  at  the  impact  face,  thereby  enabling  evidence  of  reaction 
initiation  associated  with  compaction  or  compact  deformation  processes  during  impact. 


Figure  9.  Schematic  of  sapphire  anvil  set-up  (a)  and  Imacon  high  speed  images  of  reaction  occurring  at  the 
impact  face  of  the  sapphire  anvil. 
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Figure  9  shows  a  schematic  of  the  sapphire  anvil  setup  and  images  of  an  aluminum  powder 
compact  impacting  the  sapphire  anvil.  A  frame  is  captured  every  500ns  and  the  images  show  an 
initial  impact  flash  that  is  extinguished  within  lps,  followed  by  reaction  light  that  observed  in  a 
few  localized  areas  within  the  compact.  These  images  need  to  be  correlated  with  a  profile  view 
of  the  impact  process  to  determine  at  the 
stage  of  impact  at  which  the  reaction  occurs 
during  powder  densification,  compact 
deformation,  or  post-deformation  as  the 
copper  projectile  interacts  with  the  anvil. 

Additionally,  IR  imaging  can  be  employed  in 
conjunction  with  the  sapphire  anvil  (as 
shown  in  Figure  10)  in  order  to  identify  the 
nature  of  heating  occurring  during  reaction, 
whether  localized  or  bulk.  The  analysis  of 
impact  initiation  trends  as  a  function  of 
particle  size,  porous  density,  and  mechanical 
pre-activation  can  then  be  correlated  with  in 
situ  time-resolved  analytical  methods, 
including  IR  camera  measurements  and 
UV/Vis  spectroscopy,  to  identify  the  source 
of  the  differences  in  reaction  phenomena 
between  compacts  of  aluminum  powders  of 
different  characteristics. 


Figure  10.  IR  (a-b)  and  corresponding  Imacon  camera 
(c)  images  of  aluminum  compact  impacting  sapphire 
anvil  at  267m/s.  IR  images  were  taken  during  initial 
impact  (a)  and  3ms  following  impact  (b). 


3.2  Impact-initiated  Reactions  in  Ti+Al+B  Mixtures  Under  Uniaxial  Stress  Loading 

The  objectives  of  this  task  are  to  investigate  the  reaction  behavior  of  mixtures  in  the 
Ti+Al+B  system  with  Ti+2B  mixed  with  varying  volume  fractions  of  Al  powders.  Ultimately,  we 
seek  to  understand  the  reaction  mechanisms  in  this  ternary  system  and  to  explore  the 
synergistic  effects  that  Ti+B  intermetallic  reaction  may  have  on  Al  powders  in  promoting  an 
oxidation  reaction.  Moreover,  this  work  aims  to  quantify  the  extent  of  the  role  that  meso-scale 
features  (extrinsic  properties)  play  in  ensuring  the  right  thermodynamic  conditions  are  ensured 
so  that  these  systems  may  react  in  a  deterministic  manner.  The  investigation  involved 
performing  experiments  under  uniaxial-stress  and  simulating  experimental  conditions  using 
realistic  microstructural  representations  using  CTH  hydrocode. 

Uniaxial  stress  rod-on-anvil  impact  experiments  were  performed  on  Ti+Al+B  mixtures, 
varying  in  composition  from  0%  Al  to  75%  Al  by  volume  and  mixed  with  stoichiometric  Ti+2B. 
The  powder  mixtures  were  pressed  into  1/8"  diameter  pellets  and  mounted  on  copper  rods 
which  were  impacted  against  a  steel  anvil  in  the  7.62  mm  gas  gun.  A  number  of  impact 
experiments  characterized  using  high  speed  video  and  a  framing  camera  were  conducted  to 
observe  the  reaction  event  confirmed  via  light  emission.  The  "go/no-go"  experiments  revealed 
whether  the  pressed  powder  mixture  pellet  produced  light  emission  (indicative  of  reaction  as 
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shown  in  representative  micrograph  in  Figure  11)  at  a  certain  firing  velocity  (kinetic  energy), 
which  was  used  to  determine  the  "threshold"  energy  for  reaction  initiation. 


75%  Al  -  25%  Ti+2B 
Velocity  =  454.21  m/s 

Initiation  -  5  ps  8  ps  1 1  ps 

Figure  11:  Representative  images  of  typical  light  emission  observed  in  a  rod-on-anvil  impact 
experiment,  indicative  of  possible  initiation  of  intermetallic  or  oxidation  reaction. 

It  was  found  that  the  50%AI  mixture  reacted  most  easily  under  uniaxial  stress 
conditions,  as  evidenced  by  its  lower  reaction  threshold  energy  versus  that  for  the  more 
deficient  (pure  Ti+2B)  and  Al-rich  Ti+2B  +  50  vol%  Al  mixtures.  The  enhanced  reactivity  of  this 
mixture  can  be  related  to  the  competition  between  dissipative  mechanisms  being  prevalent  in 
Al  deficient  mixtures  and  localization  being  hindered  with  higher  Al  content  mixtures.  Results  of 
the  go/no-go  reaction  experiments  are  summarized  in  Figure  12,  which  plots  impact  energy 
against  material  composition  with  pure  Al  composition  (over  90%  TMD)  added  for  comparison. 
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Figure  12:  Reaction  threshold  map  showing  the  energy  requirements  to  observe  reactivity 
under  uniaxial  stress  conditions  for  a  variety  of  stoichiometries. 


A  lower  reaction  threshold,  defined  as  the  amount  of  energy  required  to  observe  light 
emission  in  more  than  two  frames  during  impact,  is  observed  for  the  50%  Al  stoichiometry  than 
for  any  of  the  other  mixtures.  Aluminum  deficient  materials  show  a  slightly  higher  reaction 
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threshold.  The  reactions  observed  in  Al  deficient  materials  are  likely  driven  by  frictional  heating 
and  localized  hot-spot  formation  as  a  result  of  pore  collapse  during  the  compression  process. 
The  Al  rich  mixtures  rely  on  strain  localization  phenomena  to  initiate  a  reaction.  However,  for 
greater  Al  content,  the  amount  of  Ti+2B  to  satisfy  stoichiometric  requirements  in  the  mixture 
needs  to  decrease.  This  may  lead  to  any  reaction  initiation  being  quenched  before  a  sustained 
reaction  between  the  Ti+2B  can  sufficiently  ignite  the  Al  reaction.  The  lower  temperatures  in  Al- 
containing  materials  may  be  sufficient  to  enhance  diffusion-controlled  reactions  in  Al,  which 
can  trigger  the  Ti+2B  reactions.  Quantitative  stereology  has  revealed  that  a  more  intimate 
contact  between  Ti  and  B  can  be  obtained  during  mixing  in  the  presence  of  Al  content.  The 
mechanistic  effects  at  the  meso-scale,  that  affect  the  reaction  initiation  behavior  of  these 
materials  was  further  investigated  using  numerical  simulations. 

Uniaxial  stress  simulations  were  performed  on  the  Ti+2B+AI  powder  mixtures  using  real 
microstructures.  Figure  13  shows  the  simulation  domain  with  explicit  description  of  boundary 
conditions.  Results  of  simulations  performed  at 
200  m/s  reveal  substantial  agglomeration  of 
smaller  boron  particles  in  well-defined  dark 
bands  near  localization  regions  of  large  Ti  or  Al 
deformations,  as  shown  in  Figure  14.  Large 
temperatures,  in  excess  of  1200°  C  are  observed 
even  at  low  impact  velocities  in  Ti+2B  mixture 
(with  0%  Al),  largely  due  to  the  large  pore  space 
where  particles  can  mix  and  deform.  This  leads 
to  increased  friction  and  deformation  of  softer 
load-carrying  Al  and  Ti  particles  which  account 
for  the  high  temperatures.  While  local  temperatures  are  higher  in  Ti+2B,  bulk  temperature 
remains  low  which  can  account  for  higher  reaction  threshold.  Bulk  temperatures  in  Al- 
containing  structures  are  sufficient  to  cause  Al  to  melt  and  thus  possibly  lead  to  a  potential 
liquid-phase  mediated  reaction  mechanism.  Transport  phenomena  are  currently  being  studied 
using  stereological  metrics  to  assess  the  evolution  of  particle  interfaces,  agglomeration,  and  hot 
spot  regions.  The  agglomerated 
structure  is  better  visualized  in 
the  zoomed-in  view  with  the 
superimposed  temperature  field, 
as  shown  in  Figure  15.  Both,  the 
0%  Al  and  50%  Al  structures 
show  pronounced  banding  and 
agglomeration.  However,  the 
number  of  small  boron 
agglomerates  is  much  higher  for 
the  0%  Al  structure,  indicating 
that  the  presence  of  Al 
somewhat  reduces  the  boron 
agglomeration.  This  may  be  due 


Figure  14.  Half-view  of  uniaxial  stress  impact  simulation  at  200  m/s 
showing  pronounced  agglomeration  of  boron  in  bands.  Although 
local  temperatures  are  higher  with  0%  Al,  potential  melting  of  Al  in 
50%  Al  sample  can  produce  reaction. 


Sound  speed  absorbing  boundary  (half-space) 
Figure  13.  Uniaxial  stress  simulation  domain 
with  explicit  description  of  boundary  conditions. 
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to  Al  preferentially  shearing  and  filling  empty  space,  whereas  the  0%  Al  structure  depends  on 
shearing  of  Ti  and  movement  of  the  smaller  B  particles  to  fill  the  empty  spaces  during  crush-up. 
The  zoomed  view  in  figure  15 
compares  regions  near  the 
middle  of  the  structure  for 
different  times.  The  structure 
of  the  agglomerates  is  revealed 
as  colonies  of  boron  particles, 
stretched  and  in  some  cases 
pulverized  (most  prevalent  in 
the  0%  Al  structure).  The  boron 
particles  are  much  more 
equiaxed  in  the  0%  Al  structure 
than  in  the  50%  structure.  The 
Al  in  the  50%  structure  carries 
boron  along  as  it  deforms, 
whereas  the  0%  structure 
leaves  greater  open  spaces 
since  the  Ti  does  not  deform  as 
readily.  These  open  spaces  do 
not  permit  boron  deformation. 

Instead,  boron  agglomerates 
and  fractures  during  the  mixing  phase  within  these  open  spaces  (i.e.  pores).  The  bulk 
temperature  also  is  much  greater  in  the  0%  Al  structure  where  the  heating  is  based  on  the 
plastic  dissipation  in  the  Ti  and  the  interparticle  friction  between  Ti  and  B  surfaces.  These 
dissimilar  effects,  most  clearly  demonstrated  by  meso-scale  simulations  illustrate  the  effects  of 
the  reactant  configuration  on  impact-initiation  of  reactions  in  Ti+Al+B  ternary  mixtures. 


0%  Al  50%  Al 

Figure  15.  Zoomed-in  view  near  centerline  of  microstructure  showing 
bands  of  B  particles  and  higher  temperatures  in  0%  Al.  Temperatures 
for  50%  Al  are  significant  and  can  exceed  melt  temperature. 


3.3  Shock-compression  of  Ti+2B  powder  mixtures:  Experiments  and  Simulations 

The  highly  heterogeneous  starting  Ti+2B  powder  mixtures  and  the  continually  evolving 
nature  of  the  heterogeneities  can  result  in  various  reaction  pathways  which  make  tailoring  of 
the  reaction  configuration  a  challenge  to  design.  The  Ti/B  system  is  of  particular  interest 
because  of  the  extraordinarily  high  enthalpy  of  reaction  (AHrxn0  =  -320  kJ/mol)  when  forming 
the  stable  TiB2  product  phase.  This  unconventional  energetic  mixture  has  many  key  ingredients 
that  can  enable  shock-induced  reactivity  due  to  meso-scale  transport  processes  and  strain 
localization  phenomena  enhanced  by  the  very  heterogeneous  nature  of  these  mixtures.4,10 

Experiments:  The  shock  compression  response  of  1:2  molar  mixtures  of  Ti+B  was  investigated 
at  shock  pressures  up  to  5  GPa  to  probe  the  shock  adiabat  of  this  mixture  for  a  50%  starting 
density.  The  shock  compression  regime  was  probed  from  crush-up  to  full  density  and  onward  to 
assess  the  potential  onset  of  a  shock-induced  chemical  reaction  event  in  the  powder  mixture. 
PVDF  stress  gauges  measure  the  shock  speed  and  equilibrated  stress  at  the  gauge/powder 
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interface  and  the  Velocity  Interferometer  System  for  Any  Reflector  (VISAR)  was  employed  to 
compare  the  interface  velocity  at  the  powder/backer  interface  with  the  PVDF  stress  signal.  A 
series  of  two-dimensional  continuum  meso-scale  simulations  on  validated  simulated 
microstructures  were  performed  to  predict  the  shock  compression  response  and  identify  the 
meso-scale  mechanics  that  are  essential  for  reaction.  The  suitability  of  the  simulated 
microstructural  representations  was  evaluated  by  comparing  the  experimental  and  predicted 
pressure  traces.  The  measured  stress  profiles  are  compared  with  meso-scale  simulations  of  a 
highly  resolved  uniform-random  mixture  of  Ti+2B  also  at  50%  TMD.  The  measured  backer  stress 
gauge  traces  show  an  interesting  "two-wave"  rise-time  structure  stemming  from  particle 
avalanching  effect  caused  by  initial  compression  of  loosely-bound  particles  near  the  surface  of 
the  PVDF  stress  gauge,  followed  by  an  avalanche  of  fully  compressed  particles  increasing  the 
stress  to  maximum  levels.  The  simulations  are  found  to  generally  over-predict  the  measured 
stress  at  the  backer,  though  they  capture  the  salient  features  of  the  rise-time  effects  caused  by 
the  heterogeneous  powder  compaction  which  manifest  themselves  in  the  averaged  stress 
traces  as  discrete  humps  and  compare  very  favorably  with  measured  PVDF  traces. 

Uniaxial  strain  experiments  were  performed  on  mixtures  of  titanium  (Atlantic 
Equipment  Engineers  ~625  mesh  sponge)  and  boron  (Alfa-Aesar  Amorphous  94-96%,  Mg 
nominal  1%  <  5  micron  irregular/flake)  powders  first  dried  at  110°C  and  then  pre-mixed  in  1:2 
molar  ratio  in  an  inert  atmosphere.  The  powders  were  sealed  and  transferred  to  a  v-blender  for 
further  mixing  (24  h)  to  ensure  a  uniform-random  mixing  of  Ti  and  B.  Further  details  of  the 
powders  used  in  this  study  can  be  found  elsewhere.45  Instrumented  parallel-plate  impact 
experiments  on  Ti+2B  powders  were  conducted  using  the  80  mm  helium-driven  gas  gun.  A  set 
of  copper  capsule  fixtures  were  prepared  to  contain  the  powder  mixture  during  the 
experiments.  The  capsules  (Figure  16)  consisted  of  an  outer  copper  containment  ring,  copper 
driver  plate,  and  transparent  fused  silica  backer  plate  sandwiching  the  powder  which  was 
pressed  onto  a  dummy 
driver  at  50%  TMD.  The 
dummy  driver  was 
replaced  with  one 
containing  the  gauge 
package.  Aluminum 
sabots  carried  the  flyer 
plate  (e.g.,  copper, 
tungsten,  or  tungsten 
heavy  alloy  depending 
on  the  desired  impact 
stress)  and  impact 
velocity  was  monitored  by  sequential  shorting  pins  as  shown  in  the  schematic  in  Figure  17. 
PVDF  stress  gauge  packages  were  mounted  onto  the  driver  and  backer  plates  to  monitor  the 
stress  state  and  wave  arrival  times  at  the  front  and  back  of  the  powder  bed.  The  gauge 
packages  consisted  of  alternating  polytetrafluoroethylene  (PTFE)  layers  to  isolate  the  PVDF 
gauge,  and  a  sputtered  layer  of  Al  (150  nm  thick)  at  the  end  to  mitigate  the  effects  of 
pyroelectric  interference  on  the  PVDF  gauge  (Figure  16). 


Figure  16:  Capsule  assemblies  showing  driver  (Cu)  and  backer  (Fused  silica) 
with  attached  PVDF  stress  gauge  packages.  Surface  of  gauge  package  is 
sputter  coated  with  Al  (150  nm)  to  prevent  pyroelectric  interference  gauge. 
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PVDF  gauge  leads  (yellow) 


Velocity  shorting 
pins 


Figurel7.  Schematic  of  uniaxial  strain  impact  experiment  set-up.  An  Al  sabot  carries  a  copper  impactor  plate 
which  impacts  a  Cu  driver  plate  and  imparts  the  shock  into  the  powder  backed  by  a  fused  silica  backer  plate. 
PVDF  stress  gauge  packages  sandwich  the  powder  and  provide  stress  and  shock  speed  measurements. 


Figure  18  shows  an  example  of  stress  profiles  captured  by  the  input  and  propagated 
PVDF  stress  gauge  profiles.  The  PVDF  gauges  at  the  driver/powder  interface  showed  anomalous 
signals  at  the  peak  stress  in  several  experiments  possibly  due  to  gauge  failure.  Therefore,  the 
shock  input  stress  in  the  powder  mixture  was  impedance-matched  with  the  shock  speed  and 
impact  velocity  at  the  driver.  The  propagated  gauge  stress  profiles  captured  at  the  powder- 
backer  interface,  reveal 
characteristic  effects  of  wave 
dispersion  associated  with  the 
heterogeneous  nature  of  powder 
compaction.  The  effects  are  clearly 
revealed  by  the  backer 
(propagated)  stress  gauge  profiles 
shown  in  Figure  19.  The 
transparent  fused  silica  backer  also 
permitted  for  simultaneous  VISAR 
measurements  through  adequate 
preparation  of  the  fused  silica 
surface  prior  to  attaching  the 
gauge  package.  All  mating  surfaces 
were  lapped  to  optical  flatness  (<  1 
mrad  tilt). 


Figure  18.  Typical  stress  profiles  captured  by  input  and 
propagated  PVDF  stress  gauges. 
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Uniaxial-strain 
plate-impact 
experiments  provide 
the  shock  speed  and 
equilibrated  stress  at 
the  powder/backer 
interface.  The  results 
for  each  experiment 
are  summarized  in 
Table  2.  The  shock 
speed  through  the 
powder  was  calculated 
based  on  the  transit 
time  measured  at  the 
10%,  50%,  and  90% 
levels  of  the  rise  time 
to  peak  stress.  This 
provides  three 

independent 
measurements  of  shock 

speed,  though  one  can  see  considerable  scatter  due  to  the  dispersion  at  the  backer  gauge.  The 
dispersion  decreases  as  the  impact  speed  and  impactor  plate  impedance  increase.  The  stress 
waves  propagating  through  the  powder  mixture,  and  recorded  at  the  interface  with  the  backer 
plate,  result  in  significant  dispersion  at  lower  impact  velocities,  which  increases  the  transit  time 
difference  between  rise  time  levels  leading  to  increased  uncertainty  in  the  single  shock  speed 
value.  The  thickness  of  the  powder  was  measured  through  five  independent  measurements 
with  an  accuracy  of  0.005  mm.  The  greatest  contribution  to  the  uncertainty  in  thickness  arises 
from  irregularities  during  the  initial  pressing  of  the  powders  in  the  capsule.  The  overall  results 
reveal  no  clear  trend  and  no  discernable  evidence  of  a  chemical  event. 


Available  backer  stress  traces  for  successful  experiments 
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Figure  19.  Backer  stress  traces  show  a  variety  of  rise-time  behaviors 
characteristic  of  highly  heterogeneous  mixtures.  Trace  for  818  m/s  shows 
numerous  signals  which  may  be  due  to  separate  compression  phenomena 
corresponding  to  two  separate  voltage  spikes  in  raw  signals  (not  shown).  The 
higher  imapct  velocity  and  higher  impedance  flyer  experiments  show  this  same 
phenomean,  but  not  as  pronounced  due  to  lower  inherent  dispersion  in  rise 


Table  2:  Summary  of  experimentally-measured  quantities  and  impedance-matched  particle 
velocity  and  stresses  for  the  shock  compression  experiments  on  Ti+2B. 


Shot  # 

Flyer 

Vimp  (mm/ps) 

Us  (mm/ps) 

Up  (mm/ps) 

a/  (GPa) 

V/Voo 

1225 

Cu 

0.513  ±  0.35% 

1.23  ±  4.14% 

0.48  ±  7.08% 

1.06  ±  4.06% 

0.61 

1227 

Cu 

0.818  ±  1.61% 

1.76  ±  1.95% 

0.753  ±  4.54% 

2.24  ±  1.69% 

0.57 

1302 

Cu 

1.044  ±  3.66% 

2.16  ±  5.30% 

0.944  ±  3.62% 

3.63  ±  1.49% 

0.56 

1307 

Cu 

0.832  ±  0.28% 

1.72  ±  2.02% 

0.769  ±  4.45% 

2.36  ±  2.37% 

0.55 

1309 

WHA 

0.921  ±  0.52% 

1.79  ±  3.68% 

1.01  ±  3.38% 

3.83  ±  15.7% 

0.43 

1310 

W 

0.933  ±  0.52% 

2.13  ±  3.44% 

1.12  ±  3.06% 

4.1  ±  7.32% 

0.48 

1314 

W 

0.890  ±  0.66% 

1.81  ±  0.40% 

1.08  ±  3.17% 

3.84  ±  13.01% 

0.42 

1320 

Cu 

1.06  ±  0.83% 

1.86  ±  1.86% 

0.966  ±  3.54% 

3.03  ±  0.71% 

0.48 
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Simulations:  Realistic  simulated  microstructures  were  generated  by  inserting  Ti  particles 
from  a  specially-crafted  library  created  through  montage  serial  sectioning,  and  simulating 
boron  particles  through  a  polar  plot  of  a  Fourier  Sine  series.  Further  details  of  the  method  are 
described  in  a  prior  work.30  The  particles  were  inserted  into  a  simulation  domain  one  at  a  time 
with  a  strict  overlap  control.  The  size  distribution  for  each  component  replicated  the  vendor- 
supplied  information  and  a  uniform-random  distribution  of  constituents  was  generated.  A  50% 
TMD  uniform-random  distribution  of  particles  was  achieved  and  validated  through  two-point 
correlation  functions  when  compared  with  real  microstructures. 


Distance 

Figure  20:  Schematic  of  CTH  simulation  domain  for  uniaxial  strain  simulations. 


Uniaxial-strain 
impact  simulations 
were  performed  using 
the  Eulerian  hydrocode 
CTH  (Version  9.0, 

Sandia  National 

Laboratories  2009)  by 
placing  the 

microstructure  in  a  2D 
simulation  domain  with 
simulated  PVDF  gauge 
packages,  Cu  driver, 
impactor  (either  Cu  or 
W),  and  fused  silica 
backer  as  shown  in 
Figure  20.  The 
simulation  set-up  is  a 

1:4  scale  model  of  the  experimental  fixture.  The 
left  and  right  boundaries  were  set  as  periodic 
symmetry  boundaries  and  the  top  and  bottom 
boundaries  were  set  as  sound  speed  absorbing 
(half-space).  The  Ti  particles  were  modeled  with 
the  rate-independent  Steinberg-Guinan-Lund 
constitutive  model.  The  Cu  plates  were  modeled 
with  the  Zerilli-Armstrong  constitutive  model 
and  the  gauge  package  was  made  hydrodynamic. 

The  fused  silica  backer  employed  a  Johnson-Cook 
model  and  the  boron  particles  used  a  Von-Mises 
yield  surface  with  thermal  softening.  All 
constituents  except  for  boron  were  modeled 
with  the  SESAME  equation  of  state  as  provided 
by  the  built-in  SESAME  tables;  Boron  was 
modeled  with  a  Mie-Gruneisen  equation  of  state.  A  set  of  ten  tracer  points  were  placed  at  each 
PVDF  gauge  center  to  track  the  velocity,  stress,  and  temperature.  A  typical  temperature  profile 
is  shown  in  Figure  21  for  a  fully  compacted  structure. 


Temperature  (K) 
■  1400 


Figure  21.  Temperature  field  for  a  1000  m/s 
Cu/Cu  impact  on  a  50%  TMD  Ti+2B  powder. 
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The  meso-scale  microstructure-based  simulations  show  a  random  distribution  of  hot 
spots  exceeding  1400  K,  but  do  not  predict  a  large  heat  buildup  in  the  bulk  of  the  material.  This 
can  lead  to  quenching  of  any  potential  reaction.  Temperatures  remain  nominally  within  the 
700-1000K  range,  well  below  even  the  temperatures  needed  for  thermally-induced  reactivity  in 
Ti+2B.  Also,  the  diffusional  processes  within  these  temperature  ranges  are  too  sluggish  to 
promote  enhanced  reactivity.  However,  there  was  some  indication  of  shear  banding  caused  by 
excessive  deformation  in  the  Ti  particles,  potentially  leading  to  increased  temperatures  within 
the  surrounding  regions.  The  size  disparity  and  differences  in  elastic-plastic  properties  between 
Ti  and  B  particles,  causes  the  B  particles  to  agglomerate  and  confine  between  the  larger  Ti 
particles.  However,  the  high  strength  of  B  and  the  geometry  of  the  Ti  particles  do  not  permit 
large  hydrodynamic  transport  as  previously  observed.  This  can  unequivocally  be  shown  through 
the  evolution  of  surface  area  or  B/B,  Ti/Ti  contact  area  as  a  function  of  simulation  time.  Studies 
are  currently  being  conducted  by  varying  extrinsic  properties  such  as  particle  anisotropy, 
degree  of  starting  agglomeration,  and  size  distribution  using  a  2k  design  of  experiments 
framework  to  ascertain  the  sensitivity  of  the  system  to  variations  in  adjustable  parameters. 


The  recorded  backer-gauge  stress  traces  (Figure  22)  show  a  dispersed  "two-wave" 
structure,  potentially  due  to  differences  in  arrival  times  of  the  more  loosely-bound  particles 
near  the  surface  of  the  PVDF  gauge  and  the  succeeding  particle  "avalanche."  The  avalanching 
particles  move  at  a  group  velocity  impedance-matched  with  the  Cu  driver  and  slower  than  the 
initial  jagged  front.  This  initial  front  creates  an  initial  shock  compression  of  the  gauge  which 
rings  up  to  peak  pressure  for  about  10-20  ns,  subsequently  equilibrating  with  the  more  sluggish 
particles  that  follow.  This  is  confirmed  by  simulations  which  reproduce  the  same  behaviors  (c.f. 
Figures  22)  due  to  the  continual  compression  of  the  backer  by  a  fully  compacted  powder. 


Figure  22:  Simulated  (black)  and  experimentally-measured  (red)  PVDF  gauge  stress  profiles  showing  the 
correlation  for  low  and  high  velocity  experiments.  The  513  m/s  experiment  (left)  shows  re-shock  event  due  to 
re-compression  wave  emanating  from  the  fully  dense  powder  compacted  powder,  which  is  captured  both  by 
the  PVDF  gauge  and  simulations  to  a  very  high  fidelity.  In  the  case  of  the  933  m/s  experiment,  the  simulate 
trace  shows  a  much  more  pronounced  two-wave  structure  possibly  due  to  the  larger  gauge  package 
(unsealed)  employed.  However,  the  salient  features  are  captured  well.  The  rise-time  is  also  much  shorter  in 
the  simulated  trace  due  to  added  stiffness  of  the  particles  because  of  the  sliding  algorithm  imposed  bv  CTH. 
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The  simulations  also  correlate  extremely 
well  with  VISAR  traces  as  shown  in  Figure 
23,  capturing  the  peak  velocity  and  even 
release-wave  events,  providing  an  initial 
validation  for  the  simulation  scheme  and 
physical  models  employed.  The  simulated 
trace  even  illustrates  impactor  edge  release 
effects,  although  simulated  trace  shows 
earlier  release  than  the  one  captured  by 
experimental  trace.  However,  the 
simulated  peak  stresses  also  over-predict 
the  recorded  PVDF  stress  profiles.  This  may 
be  due  to  an  overly-stiff  response  from  the 
model,  or  likely  premature  failure  of  gauge. 
Further  experiments  are  needed  for 
validation. 


Particle  velocity  in  backer  \'tmpact  =  890  m/s 
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Figure  23:  Simulated  (black)  particle  velocity  profiles  for 
933  m/s  experiment  show  extremely  good  correlation 
with  non-contact  measurement  of  particle  velocity 
profile  captured  using  VISAR  (red). 


4.  Statistical  Representation,  Reconstruction,  and  Simulations  of  Meso-Scale  3-D 
Microstructures  of  Reactive  Materials 

In  the  powder-processed  ternary  Ti-B-AI  compacts,  the  relative  spatial  arrangements  of 
Ti,  B,  and  Al  particles  are  expected  to  influence  the  kinematics  of  impact  initiated  reactions  as 
well  as  deformation,  fracture,  and  fragmentation  of  the  constituent  particles  during  high  strain 
rate  processes.  However,  such  effects  have  not  been  systematically  investigated  in  the  past.  To 
understand  the  effects  of  powder  processing  techniques  on  the  microstructural  spatial 
arrangements;  to  develop  and  implement  the  techniques  to  quantify  and  mathematically 
represent  such  differences  in  the  spatial  arrangements;  to  formulate  approaches  to  incorporate 
such  variations  in  the  spatial  arrangements  in  the  simulated  realistic  multi-phase 
microstructures;  and  to  simulate  the  effects  of  such  microstructure  variations  on  the 
deformation,  fracture,  and  fragmentations  of  the  particles  during  high  strain  rate  deformation 
using  realistic  simulated  microstructures,  a  series  of  experiments  were  performed.  A  set  of  Al- 
Ti-B  specimens  having  a  range  of  compositions  were  prepared  using  two  different  powder 
mixing  routes.  In  one  set  of  specimens,  first  Ti  and  B  powders  were  pre-mixed,  and  then  in  the 
second  step,  the  pre-mixed  Ti-B  was  blended  with  Al  powder  and  compacted  to  obtain  ternary 
A-Ti-B  compact  of  different  compositions.  In  the  second  series  of  specimens,  all  three  elemental 
powders  (Al,  Ti,  and  B)  were  simultaneously  blended  and  then  compacted  to  obtain  the  ternary 
compact  of  the  same  set  of  compositions.  These  different  powder  mixing  routes  lead  to 
different  spatial  arrangements  of  the  Ti,  Al,  and  B  constituents  in  the  compacts  having  the  same 
volume  fractions  and  size  distributions  of  Al,  Ti,  and  B  phases. 
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Figure  24  illustrates  microstructures  of  two  such  specimens  having  the  same  Ti,  B,  and 
Al  volume  fractions,  size  distributions,  and  all  other  first  order  microstructural  properties.  In 
these  micrographs,  Y-axis  is  the  compaction  direction.  Both  microstructures  in  the  figure  appear 
to  be  isotropic  and  with  no  visible  differences  in  the  spatial  arrangements  of  the  three  phases. 


(a)^ 


Figure  24.  Microstructure  of  AI-4vol.%  Ti-14vol.%B  powder  compacts  (a)  Ti-B  powder  pre-mixing  followed 
by  blending  of  Al  powder  with  pre-mixed  Ti-B,  (b)  Al,  Ti,  and  B  powders  simultaneously  blended.  The  bright 
phase  is  Al,  grey  particles  are  Ti,  and  dark  phase  is  B.  Y-axis  of  the  micrographs  is  the  compaction  direction. 


However,  quantitative  measurements  of  microstructural  two-point  correlation  functions  clearly 
reveal  that  (i)  there  is  anisotropy  of  spatial  correlations  in  these  microstructures,  and  (ii)  spatial 
correlations  among  Ti,  B,  and  Al  constituents  are  indeed  different  in  the  two  microstructures. 
Figure  25  depicts  the  Ti-Ti  two  point  correlation  functions  in  the  two  powder  compacts  along 
the  compaction  direction  and  the  radial  direction  (i.e.,  perpendicular  to  compaction  direction) 
prepared  by  simultaneously  blending  the  Al,  Ti,  and  B  powders.  Interestingly,  the  two-point 
correlation  functions  in  the  two  directions  are  significantly  different. 
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Figure  25.  Ti-Ti  two-point  correlation  function  in  compacts  along  compaction  direction  (red)  and  radial 
direction  (blue).  Samples  were  prepared  by  simultaneously  blending  Al,  Ti,  and  B  powders  followed  by  uniaxial 
compaction  of  loose  powder  stack  at  room  temperature  (a)  AI-4vol.%  Ti-14vol.%B,  (b)  Al-10vol.%  Ti-10vol.%B. 

Similar  trends  have  been  observed  in  the  samples  of  other  compositions  and  processing 
conditions.  Therefore,  the  process  of  uniaxial  powder  compaction  induces  mild  anisotropy  in 
the  spatial  arrangements  of  the  particles  in  these  microstructures  although  individual  Ti  (as  well 
as  B)  particles  have  equiaxed  morphologies.  These  differences  in  the  correlation  functions  along 
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the  compaction  and  radial  directions  are  probably  caused  by  particle  rearrangements  during 
uniaxial  powder  compaction  process  and  uniaxial  deformation  of  Al  which  leads  to  plastic  flow 
of  Al  in  the  pore  spaces  of  the  powder  stack  that  brings  the  Ti  particles  closer  along  compaction 
direction.  The  extent  of  such  direction  dependent  spatial  correlations  depends  on  the 
parameters  of  the  powder  compaction  processing  such  as  the  compaction  load  and  amount  of 
compaction.  Such  direction  dependence  of  the  spatial  correlations  may  in  turn  affect  the 
performance  of  the  energetic  materials. 

Different  powder  mixing  routes  (pre-mixing  of  Ti  and  B  versus  simultaneous  mixing  of  Al, 
Ti,  and  B  powders)  also  affect  the  spatial  arrangements  and  correlations  represented  by  the 
correlation  functions.  Figure  26  depicts  normalized1  Ti-B  two-point  correlation  functions  for  the 
AI-4vol.%  Ti-14vol.%B  compacts  along  the  compaction  direction  at  short  range  and 
intermediate  range  and  Figure  27  depicts  the  corresponding  behavior  along  the  radial  direction. 


1  Normalized  Ti-B  two-point  correlation  function  is  obtained  by  dividing  the  absolute  two-point  correlation  function 
by  the  product  of  the  volume  fraction  of  Ti  and  B,  fTi  and  fB.  The  normalized  Ti-B  correlation  function  approaches 
one  as  the  distance  r  goes  to  infinity 
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Note  that  the  two-point  correlation  functions  nicely  capture  the  differences  in  the  Ti-B 
spatial  correlations  resulting  from  the  two  different  powder  mixing  routes:  in  the  compaction 
as  well  as  radial  direction  the  normalized  two-point  function  probability  for  the  "pre-mixed" 
compact  (red  curve)  is  consistently  lower  than  that  for  the  compact  where  the  Al,  Ti,  and  B 
powders  are  simultaneously  mixed  (blue  curve)  up  to  the  distances  of  20  pm  indicating 
enhanced  repulsion  between  Ti  and  B  at  short  distances.  This  is  likely  to  be  due  to  larger  extent 
of  agglomeration  of  the  fine  B  particles  in  the  "pre-mixed"  compact.  The  data  on  B-B,  Ti-Ti,  and 
Al-B  two-point  correlation  functions  and  lineal  path  probability  distributions  are  consistent  with 
these  trends.  Such  data  sets  are  crucial  for  development  and  validation  of  the  techniques  for 
simulations  of  realistic  multi-phase  3D  microstructures  where  the  simulations  can  mimic  subtle 
variations  in  the  spatial  correlations  of  the  constituents  arising  from  the  variations  in  the 
processing  conditions  such  as  powder  mixing  routes. 

As  noted  in  the  above  results,  for  given  constituent  elemental  powders,  the  relative 
spatial  arrangements  of  Ti,  B,  and  Al  particles  are  expected  to  influence  the  kinematics  of 
impact  initiated  reactions  as  well  as  deformation,  fracture,  and  fragmentation  of  the 
constituent  particles  during  high-strain-rate  processes.  The  size  and  shape  characteristics  of  the 
powder  also  influence  the  said  processes.  We  have  also  developed  a  complete  set  of  tools  to 
simulate  realistic  microstructures  in  three  dimensions  for  the  full  range  of  processing  conditions 
in  use  in  impact  experiments  (i.e.  varying  compositions,  varying  anisotropy  as  determined  by 
the  extent  of  densification,  and  for  two  different  sets  of  powders).  These  simulated 
microstructures  have  the  same  statistics  (within  measurement  error)  as  the  real  specimens, 
thus  suggesting  that  the  simulations  produce  reliable  approximations  of  the  real 
microstructures  for  use  in  impact  simulations.  In  addition,  our  microstructure-simulation 
toolbox  can  simulate  hypothetical  microstructures  that  differ  from  the  realistic  microstructures 
in  one  or  more  parameters  that  closely  correspond  to  intelligible  shape  characteristics  of  the 
microstructure  or  spatial  arrangements  in  it. 
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4.1  Characteristics  of  real  microstructures 


After  careful  examination  of  all  available 
real  microstructures  the  following  three 
salient  conclusions  emerge:  (i)  all 
microstructures  exhibit  significant  long- 
range  heterogeneities;  (ii)  there  are  no 
statistically  significant  differences  between 
the  microstructures  of  premixed  powder- 
mixture  compacts  (i.e.  when  Ti  and  B  are 
mixed  together  first  and  then  blended  with 
Al)  and  non-premixed  compacts  (i.e.  when 
all  powders  are  mixed  together  at  the  same 
time);  and  (iii)  all  microstructures  show 
marked  anisotropy,  such  that  there  is  a 
significant  difference  between  the  direction 
parallel  and  perpendicular  to  loading.  To 
illustrate  the  point  about  heterogeneities, 
Figure  28  shows  two  examples  each  from  an 
old-powder  specimen2  and  a  new-powder 
specimen,  illustrating  the  one-standard- 
deviation  range  over  which  the  two-point 
correlation  function  lies  when  calculated 
over  tens  of  different  small  microstructure 
windows  taken  out  of  a  large  montage.  Only 
such  a  calculation  can  reveal  long-range 
heterogeneities  in  the  microstructure.  Since 
the  correlation  function  at  any  given  r  is 
seen  to  follow  a  normal  distribution,  there  is 
a  roughly  a  68%  chance  that  a  randomly 
chosen  window  of  the  size  indicated  out  of 
the  montage  will  have  a  correlation  function 
lying  within  the  band  shown. 

Simulated  microstructures  (which 
typically  have  similar  window  sizes)  can  only 
be  expected  to  match  the  real 
microstructure  in  two-point  correlations  to 
about  the  same  error.  Fig.  29  plots  the 
correlation  functions  along  two  directions 
(parallel  to,  and  perpendicular  to,  loading) 


2  Throughout,  “old  powders”  refers  to  crystalline  titanium  (Cerac  325 -mesh)  and  boron  particles  (Alfa-Aesar),  and 
“new  powders”  refers  to  sponge  (crystalline)  titanium  (~  20  pm)  and  amorphous  boron  (<  5  pm).  The  aluminum 
used  for  the  specimens  was  Valimet  H2  (<  6  pm). 
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Figure  28.  One-sigma  error  bounds  on  (a)  P(Ti  Ti),  (b) 
P (B-B)?  for  sample  made  with  old  powders,  and  on  (c) 
P (B-B)  (d)  P (Ti— Ti)  for  sample  made  with  new  powders. 
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on  the  same  plot  of  a  microstructure  of  (a)  old- 
powder  compact  and  (b)  new-powder 
compact.  The  parallel-direction  correlation 
function  consistently  has  lower  slope  at  0,  as 
would  be  expected  if  particles  reoriented 
during  compression,  in  every  microstructure 
examined. 


4.2  The  Simulation  Method 


Simulations  use  real  particle 
morphologies  extracted  from  reconstructions 
of  volumes  of  real  binary  microstructures  (Al— 

Ti  and  Al-B),  with  the  exception  of  new  boron 
particles  that  proved  too  small  to  be  accurately 
segmented  even  at  the  highest  resolutions  in 
optical  microscopy.  Hence,  the  boron  particle 
shapes  are  simulated  based  on  observed  and 
known  size  and  shape  characteristics.  The 
same  problem  prevented  us  from  obtaining 
accurate  segmentations  of  ternary  images  for 
characterization.  Fig.  30  shows  some  examples 
of  Ti  and  B  particles  from  this  new  library.  The 
new  library  supplements  the  old  library  of  particles  of  old  powders  of  Ti  and  B. 
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Figure  29  Anisotropy  in  (a)  Old  powder  compact 
and  (b)  New  powder  Compact.  Blue  curve  is  parallel 
to  loading  and  red  curve  is  perpendicular  to  loading. 
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Figure  30:  (a)  New  Ti  particle  morphologies,  extracted  from  reconstructed  micro  structural  volumes.  Sizes  shown 
are  diameters  of  volume-equivalent  spheres,  (b)  New  B  particle  morphologies,  simulated  based  on  known  and 
observed  size  and  shape  characteristics,  shown  at  resolution  equal  to  that  in  the  Ti  particles  in  (a).  Sizes  shown 
are  diameters  of  volume-equivalent  spheres. 


The  simulation  method  developed  is  a  Monte  Carlo  method,  running  a  Metropolis 
simulation  within  a  simulated-annealing  framework.  It  starts  out  with  loading  the  appropriate 
libraries  of  particles  of  the  right  dimension  (since  simulations  can  be  for  2-D  or  3-D 
microstructures),  asks  for  the  target  volume  fractions  and  particle-selection  criteria,  if  any 
(based  on  size  and/or  shape).  After  selecting  the  set  of  particles  (a  random  set  that  satisfies  all 
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the  specifications),  it  asks  for  the  target  distance  metric,  which  controls  clustering.  This  metric  is 
ignored  for  real-microstructure  simulation  (more  on  this  below).  Then  it  asks  if  the  particles 
need  to  be  reoriented  to  create  a  certain  level  of  anisotropy.  We  found  that  the  anisotropic  real 
microstructures  could  be  simulated  simply  by  rotating  the  set  of  particles  by  a  set  amount 
calculated  as  a  linear  function  of  porosity.  The  rationale  for  linking  anisotropy  to  porosity  is  that 
a  greater  porosity  implies  less  particle  rotation  during  compaction,  since  the  height  of  the  loose 
powder  mixture  is  reduced  by  a  smaller  amount.  From  measurements,  we  know  that  full 
densification  (achieving  roughly  95%  of  theoretical  density)  requires  about  40%  reduction  in 
height  of  the  loose  powder  during  compaction.  We  surmise  that  the  degree  of  particle  rotation 
is  a  linear  function  of  the  fractional  reduction  in  height.  Since  no  microstructures  can  be 
obtained  for  highly  porous  samples  (due  to  lack  of  structural  integrity  for  metallographic 
observation),  this  assumption  cannot  be  experimentally  tested.  The  Metropolis  simulation 
places  particles  on  a  grid  at  the  start,  and  through  random  translations  from  that  initial 
configuration,  the  microstructure  evolves  to  meet  user  specifications.  The  translations  initially 
are  purely  random  (i.e.  there  is  a  high  "temperature"  in  the  simulated  anneal),  but  then 
become  increasingly  the  kind  that  lower  the  "energy"  (a  function  of  the  particle  positions  that 
assigns  lower  values  to  configurations  closer  to  user  specification)  by  tending  to  reject 
"energy"-raising  translations.  The  energy  function  for  simulating  real  microstructures  is  one 
that  penalizes  overlaps;  for  simulating  clustered  microstructures,  it  also  penalizes  deviation 
from  the  desired  distance  metric  (more  on  this  later).  The  simulation  stops  when  "energy"  has 
reached  a  low  plateau.  Fig.  31  shows  two  example  of  3-D  microstructures— one  old-powder  and 
one  new-powder  compact.  Fig.  32  shows  a  comparison  of  two-point  correlation  functions 
between  simulated 
and  real 

microstructures  for 
one  composition  each 
of  old-powder  and 
new-powder 
specimens.  The  good 
match  suggests  that 
the  simulated 

microstructures  are 
good  approximations 
to  the  real 
microstructures. 


Figure  31:  Realistic  microstructures  containing  (a)  old  powders:  (60  pm)2  *  34 
pm  (18%Ti-12%B);  and  (b)  new  powders:  (60  pm)3  (10%Ti-10%B). 
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Figure  32:  (a)  Two-point  correlation  functions  (B-B)  of  real  (blue  band)  and  simulated  (red  line) 
microstructures  for  an  old-powder  compact  (18%Ti-12%B),  parallel  to  (left)  and  perpendicular  to  (right) 
loading;  and  (b)  Two-point  correlation  functions  (Ti-Ti)  of  real  (blue  band)  and  simulated  (red  line) 
microstructures  for  an  old-powder  compact  (18%Ti-12%B),  parallel  to  and  perpendicular  to  loading,  for  the 
same  sample  as  in  (a). 


4.3  Parameterization  of  Shape 

Invariant  Moments:  Since  material  behavior  is  clearly  sensitive  to  the  shapes  of 
constituents,  it  is  very  helpful  to  control  parametrically  the  shapes  of  powders  in  the  specimen. 
In  the  realistic  simulated  microstructures,  the  shapes  must  be  randomly  chosen  from  the 
particle  library.  But  in  hypothetical  microstructures,  it  is  helpful  to  select  them  by  some  shape- 
sensitive  criteria.  Toward  this  end,  we  added  shape-selection  capability  to  the  simulation 
toolbox  that  uses  seven  invariant  moments  to  characterize  a  particle's  shape  (calculated  from 
the  set  of  central  moments  of  the  particle  image  up  to  order  three).  A  given  shape  (such  as  a 
circle  or  an  ellipse  with  a  certain  eccentricity)  always  has  the  same  moments,  regardless  of  the 
size  or  orientation.  Though  the  numbers  are  themselves  not  easily  interpretable,  each  set  of 
moments  (thus  each  particle  shape)  is  a  point  in  7-dimensional  space,  and  the  set  of  these 
points  can  be  subjected  to  standard  clustering  analysis  to  classify  the  available  shapes  into  n 
classes,  where  the  user  chooses  n.  In  the  simulation,  rules  can  be  specified  as  to  which  shapes 
and  in  what  proportion  should  be  used  in  the  simulation.  Fig.  33  shows  three  Ti  particles  each 
from  three  of  the  shape  classes  obtained  by  setting  n  to  5. 
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Circularization/Spheroidization  of  Shapes: 

We  also  added  the  ability  to  change  existing 
shapes  to  make  them  some  percentage 
circular  (2-D  particles)  or  spherical  (3-D 
particles),  based  on  image  processing  relying 
on  distance  transforms.  The  original  shapes 
are  considered  0%  circular/spherical,  and  a 
circle/sphere  of  the  same  area/volume  as  the 
particle  represents  100%  on  the  scale.  This 
provides  a  way  to  test,  say,  the  influence  of 
flakiness  of  shape  on  deformation  and 
reaction  behavior,  since  flake  particles  would 
have  distinctly  different  invariant  moments. 
Going  from  0%  to  100%  circular,  the 
perimeter  (or  surface  area  in  3-D)  decreases. 
Fig.  34  shows  some  points  on  the  scale  for  a 
randomly  chosen  2-D  titanium  particle. 


Class  1  Class  3  Class  5 


Figure  33:  (a)  Three  particles  each  from  three  of  the 
five  shape  classes  of  Ti  particles  (new  powder), 
obtained  by  clustering  analysis  of  invariant  moments. 
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Figure.  34.  A  randomly  chosen  Ti  particle,  circularized  to  varying  degrees,  starting  with  original  particle  (leftmost) 
to  90%  circular  (extreme  right)  in  increments  of  10%.  The  area  is  held  constant. 


4.4  Parameterization  of  Clustering 

For  a  microstructure  consisting  of  various  shapes  and  sizes,  an  nth-nearest-neighbor 
distribution  of  centroids  does  not  adequately  describe  actual  surface-surface  distances 
between  particles.  We  therefore  formulated  two  metrics:  minimum  A-to-B  and  minimum  B-to- 
A  distance  distributions  for  two  particle  phases  A  and  B.  The  >4-to-B  distribution  is  the 
probability  density  function  of 
random  variable  representing  the 
distance  to  nearest  surface  point 
on  B  from  a  randomly  chosen 
surface  point  on  A;  likewise  for  B- 
to-A.  Moments  of  these 
distributions  can  then  be  used  to 
characterize  clustering  or  as 
targets  in  simulation  by 
incorporating  them  into  the 
"energy"  function.  Fig.  35  shows 
two  microstructures  simulated  in 
this  manner,  with  the  difference  in 


Figure  35.  Two  simulated  microstructures  with  different  target 
distance  metrics,  313.5  pm2  in  area.  White  =  Ti;  gray  =  boron. 
(Left)  Mean  Ti-to-B  distance  (left)  =  7.6  pm,  Mean  B-to-Ti 
distance  =  7.8  pm.  (Right)  Mean  Ti-to-B  distance  (left)  =11.0  pm, 
Mean  B-to-Ti  distance  =  10.2  pm. 
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clustering  apparent.  All  parameterizations  are  now  available  in  3-D  and  2-D,  although  it  is  easier 
to  illustrate  them  in  two  dimensions.  They  may  be  used  in  statistical  design  of  experiments 
(using  impact  simulations)  to  disentangle  the  effects  of  various  microstructural  attributes. 

4.5  Design  of  Experiments  (DOE)  for  Delineating  Effects  of  Microstructural  Parameters 

The  parametric  investigation  of  microstructural  characteristics  on  impact  behavior  was 
extended  through  a  design  of  experiments  effort  to  delineate  the  effects  of  microstructural 
parameters.  We  have  simulated  a  set  of  realistic  and  hypothetical  microstructures  for 
parametric  investigation,  from  which  examples  follow.  In  Fig  36(a)  is  a  (94  pm)2  realistic 
simulated  microstructure  with  crystalline  Ti  and  amorphous  B  particles,  10%  by  volume  of  each. 
(Titanium  is  light  gray,  boron  dark  gray,  and  aluminum  white.)  Fig  36(b-c)  are  hypothetical 
microstructures,  as  described  below,  of  the  same  composition.  For  a  microstructure  consisting 
of  various  shapes  and  various  sizes,  an  nth-nearest-neighbor  distribution  of  centroids  does  not 
adequately  describe  actual  surface-surface  distances  between  particles.  We  therefore 
formulated  two  metrics:  the  minimum  A-to-B  and  the  minimum  B-to-A  distance  distributions 
for  two  particle  phases  A  and  B.  The  A-to-B  distribution  is  the  probability  density  function  of  the 
random  variable  representing  the  distance  to  the  nearest  surface  point  on  B  from  a  randomly 
chosen  surface  point  on  A;  likewise  for  B-to-A.  The  realistic  microstructure  in  Fig  36(a)  has 
mean  Ti-to-B  distance  1.8  pm  and  mean  B-to-Ti  distance  5.8  pm.  In  Fig  36(b)  is  a  hypothetical 
microstructure  with  the  mean  Ti-to-B  and  B-to-Ti  distances  3.2  pm  and  8.3  pm  by  design.  We 
also  used  a  morphological  parameter  to  change:  circularity.  Our  simulation  toolbox  can  change 
existing  shapes  to  make  them  some  percentage  of  the  2-D  particles  circular,  based  on  image 
processing  relying  on  distance  transforms.  The  original  shapes  are  considered  0%  circular,  and  a 
circle  of  the  same  area  has  circularity  100%.  This  provides  a  way  to  test,  say,  the  influence  of 
flakiness  of  shape  on  impact  behavior.  Going  from  0%  to  100%  circular,  the  perimeter  (surface 
area  in  three  dimensions)  decreases.  Fig.  36(c)  is  exactly  the  realistic  microstructure  of  Fig. 
36(a),  except  that  all  particles  are  now  50%  circular. 
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Numerical  Proof  of  Anisotropy  of  Compact  Microstructures:  We  rigorously  analyzed 
the  real  2-D  compact  microstructures  for  quantification  of  its  features.  Quantitative  proof  of 
anisotropy  between  the  loading  direction  and  the  perpendicular  direction,  based  on  f-tests  on 
two-point  correlations  in  perpendicular  directions,  is  now  available.  Fig  37  shows  (a)  the 


IL(0)/p. m  1  IL(0)/  urn  1  Il(0)/u m-1 


Figure  37.  (a)  Micrographs  of  Al-Ti-B  green  powder  compacts  of  (1)  4.6%  Ti-15%  B-AI,  (2)  12.7%  Ti-10%  B-AI,  and 
(3)  17.7%  Ti-5%  B-AI.  Ti  is  a  light-gray,  boron  is  dark  and  very  small  particulate  phase,  and  aluminum  is  white 
matrix;  (b)  Ti-Ti  two-point  correlation  functions  PTi-Tl  parallel  to  loading  (vertical  in  micrographs)  and 
perpendicular  (horizontal  in  micrographs)  to  loading;(c)  The  p-value  as  a  function  of  r  for  the  t- test  for  equality  of 
P-n-Ti(r)  along  and  perpendicular  to  loading,  based  on  sampled  (470  pm)2  windows;  (d)  Intersections  divided  by 
length  for  Ti  surfaces  as  function  of  theta,  /L(0),  for  the  three  samples. 
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microstructures  of  three  specimens  made  with  the  same  elemental  powders,  (b)  the  two-point 
correlation  functions  in  the  two  directions,  (c)  the  t- test  as  a  function  of  distance,  and  (d)  the 
intersections  per  length  as  a  function  of  angle  made  with  the  perpendicular  to  the  loading 
direction.  Since  the  two-point  correlation  functions  are  mathematically  constrained  to  be  equal 
at  0,  the  t- test  fails  to  reject  the  hypothesis  that  the  two  functions  are  equal  there.  But  the  p- 
value  quickly  drops  to  nearly  0  (i.e.,  the  functions  are  unequal)  with  increasing  r.  This  is  exactly 
what  is  expected  if  there  is  anisotropy  of  Ti  surfaces,  because  at  low  r,  the  two-point  correlation 
function  is  dominated  by  the  morphology  of  the  microstructure  and  at  large  r  by  spatial 
features.  We  have  submitted  a  manuscript  for  publication  in  Scripta  Materialia  summarizing 
these  findings  and  reporting  on  possible  implications  for  processing  of  powders. 


5.  Non-equilibrium  Molecular  Dynamics  Simulations  of  Shocked  Ni/Al  Nanolaminates 
5.1  Simulation  of  Ni+Al  Nanolaminates 

Non-equilibrium  molecular  dynamics  simulations  were  conducted  on  Ni/Al  laminates  in 
order  to  investigate,  at  the  atomic  level,  a  system  that  has  been  well-studied  experimentally  by 
Specht  and  Thadhani.46  Simulations  were  performed  using  the  SPaSM  (Scalable  Parallel  Short- 
range  Molecular  Dynamics)  software  developed  at  Los  Alamos,  using  a  time  step  of  0.2  fs. 
Interatomic  forces  were  computed  using  the  embedded  atom  method  (EAM)  potential 
developed  by  Purja  Pun  and  Mishin  in  2009. 47  The  shock  (longitudinal)  dimension  is  z,  and 
periodic  boundary  conditions  were  employed  in  the  x  and  y  dimensions.  The  model  systems 
were  constructed  by  building  both  the  Ni  and  Al  fee  layers  with  the  111  crystallographic 
direction  parallel  to  the  z  axis.  After  minimizing  the  potential  energy,  the  system  was 
equilibrated  at  300  K  for  0.6  ps.  The  NEMD  simulation  was  performed  with  an  atomistic  piston 
located  at  z  =  0,  with  a  thickness  approximately  equal  to  the  cutoff  distance  (6.3  A)  of  the 
potential.  The  initial  velocities  of  the  remaining  atoms  are  the  equilibrated  velocities  minus  the 
piston  velocity  (i.e.,  the  piston  is  stationary  and  all  atoms  initially  move  in  the  -z  direction).  The 
piston  velocities  ranged  from  0.5  to  1.0  km/s. 

The  systems  studied  consisted  of  a  1:1  stoichiometric  ratio  of  Ni  and  Al,  with  alternating 
layers  of  the  two  metals.  The  interfaces  of  the  layers  were  planar,  and  we  varied  the  angle  (0° 
and  90°)  between  the  interfaces  and  the  shock  direction.  The  perpendicular  (_L)  simulations  are 
similar  to  simulations  performed  by  other  researchers,48  while  simulations  with  the  shock 
running  parallel  (|  |)  to  the  interfaces  have  not  been  conducted  previously.  For  the  _L 
simulations,  we  used  two  system  configurations  each  containing  ~13.4  million  atoms.  The  first 
configuration  consisted  of  3  bilayers,  with  the  Ni  layer  313  nm  thick  and  the  Al  layer  496  nm 
thick.  The  second  configuration  consisted  of  6  bilayers,  each  one  half  as  thick  as  in  the  first 
configuration.  In  the  1 1  case,  two  system  sizes  were  used  in  order  to  establish  the  size 
dependence  of  the  results.  The  smaller  1 1  system  contained  61.6  million  atoms,  and  the  larger 
1 1  system  contained  255.7  million  atoms.  Both  configurations  contained  a  single  bilayer,  which, 


36 


in  the  larger  system,  consisted  of  a  393-nm  thick  Ni  layer  and  a  620-nm  thick  Al  layer.  More 
information  on  the  initial  system  dimensions  is  given  in  Table  III. 


Table  III.  System  dimensions. 


Interf 

layer 

Ni 

Al 

thickn 

thickn 

length 

ace 

thickn 

"bHa 

N.  t 

interf 

layer 

layer 

ess  x 

ess  y 

z 

#  Ni 

#  Al 

direct 

ess 

yers 

aces 

thickn 

thickn 

dimen 

dimen 

dimen 

atoms 

atoms 

ion 

ratio3 

ess 

ess 

sion 

sion 

sion 

nm 

Nm 

nm 

Nm 

nm 

1 

1.0 

2 

3 

78.22 

a 

123.9 

2 a 

6.67 

11.56 

404.28 

1,119,7 

44 

1,117,0 

71 

1 

2.0 

6 

n 

156.4 

4a 

247.8 

4a 

6.67 

11.56 

2425.6 

7 

6,718,4 

64 

6,702,4 

30 

1 

4.0 

3 

5 

312.8 

8a 

495.6 

7a 

6.67 

11.56 

2425.6 

7 

6,718,4 

64 

6,702,4 

30 

|| 

0.6 

1 

2 

47.40 

b 

74.73 

b 

122.14 

11.81 

595.15 

30,164, 

400 

31,449, 

600 

ii 

5.0 

1 

2 

393.4 

3  b 

620.2 

9  b 

1013.7 

3 

11.81 

297.58 

125,182 

,260 

130,515 

,840 

a  example:  the  larger  of  the  1 1  systems  has  layers  5  times  as  thick  as  the  smallest  _L  system 
b  z  dimension 
cx  dimension 


For  some  small  systems,  simulations  were  continued  after  the  shock  state  was  reached. 
This  is  often  referred  to  as  a  'cook-off'  simulation  and  is  carried  out  at  fixed  volume  and  total 
energy  (micro-canonical  ensemble).  The  initial  conditions  for  the  cook-off  runs  were  taken  from 
the  NEMD  shock  trajectories  at  the  point  in  time  at  which  the  shock  wave  reaches  the  free 
surface  of  the  material.  The  precise  trajectory  point  chosen  for  the  initial  conditions  was  the 
time  slice  at  which  the  longitudinal  acceleration  for  the  final  0.66  nm  of  material  was  closest  to 
zero.  The  material  comprising  the  piston  in  the  shock  run  was  kept  frozen  in  the  cook-off 
simulation,  and  an  additional  6.3  A  of  material  was  frozen  at  the  other  (free)  surface.  Hence, 
the  system  modeled  in  the  cook-off  simulation  is  the  shocked  state  of  the  original  system  with 
6.3  A  of  material  clamped  at  each  longitudinal  end.  Our  procedure  follows  closely  the  method 
developed  by  Sewell  and  coworkers.49 

Results  and  Discussion:  A  summary  of  our  NEMD  simulation  results  is  given  in  Table  IV. 
The  final  time  is  the  time  at  which  the  fully  shocked  state  is  reached  (when  the  shock  wave 
reaches  the  free  surface  opposite  to  the  piston).  Note  that  the  final  pressure  is  also  the  peak 
pressure.  The  final  two  columns  in  Table  IV  give  a  measure  of  the  degree  of  alloying  (mixing  of 
Ni  and  Al)  that  occurs.  Ni/AI  mixing  is  quite  small  in  our  simulations  due  to  the  short  timescales 
and  the  fact  that  the  initial  model  system  is  composed  of  perfect  crystals  and  planar 
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intermetallic  interfaces  with  relatively  low  potential  energies.  More  details  on  Ni/Al  mixing  will 
be  presented  later  in  this  section. 


Table  II.  Results  from  NEMD  simulations  on  Ni/Al  nano-laminates. 
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a  The  final  pressure  is  also  the  peak  pressure.  b  Overall  fraction  of  atoms  that  have  at  least  1 
neighbor  of  the  opposite  metal.  The  value  is  the  fraction  of  the  total  Ni  or  Al  atoms  in  the 
system  divided  by  the  total  interfacial  area. 

*  configuration  is  piston/Al/Ni/Al/Ni/Al/Ni  instead  of  piston/Ni/Al/Ni/Al/Ni/Al 


Figure  38  shows  the  temperature 
plotted  as  a  function  of  time  for  the  _L 
simulations.  Separate  simulations  were 
conducted  for  piston/Ni/Al/...  and 
piston/Al/Ni/...  configurations,  and  the 
mean  temperature  for  these  configurations 
is  shown  in  Figure  38  as  well.  Due  to  the 
large  difference  in  the  compressibility  of  Ni 
and  Al,  the  simulation  results  for  these  two 
configurations  show  striking  differences 
that  will  be  discussed  in  more  detail  below. 
In  general,  the  peak  temperature  as  well  as 
the  final  temperature  of  the  fully  shocked 
state  is  higher  when  the  piston  is  adjacent 
to  a  Ni  layer.  To  extend  the  timescale  of  the 
simulation,  we  clamped  smaller  systems  (2 
bilayers,  total  length  in  z  dimension  was  400 


Fig.  38.  Temperature  profiles  plotted  as  a  function  of 
time  for  4  systems  in  the  _L  configuration.  Two  systems 
had  3  bilayers,  1  with  Ni  adjacent  to  the  piston  (dark 
red  dashed  line)  and  the  other  with  Al  adjacent  to  the 
piston  (blue  dashed  line).  The  average  of  these  two 
temperatures  is  shown  in  the  thick  black  dashed  line. 
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nm)  obtained  from  NEMD  in  their  shocked  states  and  continued  the  simulations  using 
equilibrium  MD.  Results  are  shown  in  Figure  39.  We  did  not  observe  a  significant  increase  in 
reactivity  up  to  a  nanosecond  of  simulation  time.  The  temperatures  from  _L  and  1 1  simulations 
are  plotted  together  in  Figure  40  as  a  function  of  a  dimensionless  quantity  equal  to  the  elapsed 
time  multiplied  by  the  strain  rate  for  the  simulation.  This  allows  a  direct  comparison  of 
simulations  with  variable  longitudinal  dimensions.  It  can  be  seen  that  the  orange,  green  and 
blue  lines  have  similar  temperature  plots,  while  the  gray  stays  hot  longer.  The  blue  line 
corresponds  to  a  system  with  much  thinner  Ni  and  Al  layers  than  the  other  systems.  The  plot 
shows  that  the  timescales  for  heating  as  well  as  the  peak  temperatures  are  similar  for  the  _L  and 
1 1  cases  if  the  layer  thicknesses  are  also  similar. 


Time  /  ps 


Fig.  39.  Temperature  plotted  as  function  of  simula¬ 
tion  time  for  three  systems  containing  2  bilayers  in 
-L  configuration.  Red  and  blue  lines  correspond  to 
systems  both  with  piston/Ni/AI/Ni/AI  configurations 
but  different  particle  velocities  (0.5  and  1.0  km/s, 
respectively).  Green  line  represents  a  system  with 
piston/AI/Ni/AI/Ni  configuration.  Solid  lines  are  for 
systems  having  6  bilayers  half  as  thick,  with  solid 
gray  line  being  average  of  the  two  temperatures. 


Fig.  40.  Temperature  plotted  as  a  function  of  elapsed 
time  multiplied  by  the  strain  rate  for  the  system  (the 
abscissa  is  therefore  unitless).  The  green  (orange)  line 
is  the  average  temperature  for  two  simulations,  both 
with  3  (6)  bilayers  but  one  with  Ni  adjacent  to  the 
piston  and  one  with  Al  adjacent  to  the  piston.  The 
blue  and  gray  lines  correspond  to  simulations  of 
systems  in  the  1 1  configuration.  See  Table  I  for  more 
information. 


Figure  41a  shows  shear  stress,  temperature  and  longitudinal  particle  velocity  profiles  as 
a  function  of  longitudinal  dimension  for  the  _L  simulations,  respectively.  The  configuration  of 
the  layers  is  piston/Ni/Al/Ni/Al/Ni/Al.  The  local  quantities  were  computed  in  0.66  nm  bins 
along  the  z  dimension.  In  the  first  snapshot  (t  =  25  ps),  the  wave  is  about  halfway  through  the 
first  Ni  layer  and  the  shear  stress  is  20  GPa  and  uniform  as  a  function  of  distance  behind  the 
shock  front.  The  shear  stress  in  the  first  Ni  layer  is  partially  relieved  by  a  reflection  wave  which 
begins  to  propagate  back  toward  the  piston  when  the  shock  wave  reaches  the  first  Ni/Al 
interface.  After  100  ps  the  wave  reaches  the  second  interface  (located  at  ~700  nm).  Thereafter, 
as  the  longitudinal  particle  velocity  profiles  of  Figure  37  show,  reflection  waves  become  more 
important  as  a  significant  fraction  of  the  Al  atoms  in  bilayer  1  are  moving  in  the  opposite 
direction  of  the  shock  front  at  t  =  125  and  150  ps.  The  temperature  increase  in  this  layer 
follows  the  decrease  in  stress,  and  at  the  end  of  simulation  (not  shown)  the  peak  temperature 
in  the  Al  layer  of  bilayer  1  is  almost  800  K.  In  contrast,  in  the  Al  of  bilayer  2,  the  shear  stress  is 
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not  relieved  to  the  same  degree  on  the  timescale  of  the  simulation  and  the  temperature  does 
not  climb  above  375  K.  Figure  42b  shows  the  same  simulation  results  but  for  a  system  with  a 
piston/Al/Ni/Al/Ni/Al/Ni  configuration  of  layers.  In  this  simulation,  a  layer  of  Al  is  compressed 
earliest,  and  it  appears  that  due  to  the  significantly  higher  compressibility  of  Al  compared  to  Ni, 
the  overall  response  of  the  system  to  the  shock  is  much  less  dramatic  than  in  the  simulation 
described  in  Figure  42a.  The  most  striking  difference  between  the  simulation  results  is  in  the 
temperature  profiles  at  300  ps.  While  the  simulation  with  Ni  closest  to  the  piston  (Fig.  42a) 
yields  a  peak  Al  temperature  of  800  K,  the  same  value  for  the  system  with  Al  closest  to  the 
piston  (Fig.  42b)  is  less  400  K. 
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Snapshots  of  the  local  shear  stress  and  temperature  for  a  1 1  simulation  are  shown  in  the 
x-z  plane  in  Figure  43.  As  with  all  MD  simulations  reported  here,  the  simulation  was  conducted 
with  periodic  boundary  conditions  in  the  transverse  dimensions.  Heat  flares  from  the  interface 
into  the  Al,  then  diffuses  and  eventually  temperature  becomes  more  isotropic  in  the  individual 
layers.  This  heating  occurs  with  a  very  small  of  amount  of  Ni/Al  mixing. 
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Figure  43.  Shear  stress  and  temperature  profiles  along  an  NEMD  trajectory  (thickness  in  x  dimension  is  122  nm)  as 
a  function  of  longitudinal  and  transverse  dimensions.  The  piston  is  on  the  left  at  z  =  0,  and  the  piston  velocity  was 
0.5  km/s.  The  snapshots  were  taken  at  a)  0  ps,  b)  16.8  ps,  c)  33.6  ps,  d)  50.4  ps,  e)  67.2  ps,  and  f)  84.0  ps. 


The  temperature  profiles  for  a  simulation  with 
transverse  dimension  ~1  micron  in  thick  (Fig.  44) 
shows  that  heat  'jets'  from  the  interface  laterally 
into  the  Ni  layer.  High  shear  stress  begins  when 
adjacent  jets  interact,  which  suggests  a 
dependence  of  reaction  rate  on  bilayer  thickness. 

Alloying  was  quantified  from  our  results 
by  monitoring  the  identity  of  the  nearest 
neighbors  of  each  atom  in  the  system.  The  final 
fraction  of  atoms  that  mixes,  as  well  as  the 
percent  increase  from  the  start  of  the  simulation 
is  listed  for  each  trajectory  in  Table  III.  Figure  45 
shows  the  time  dependence  of  another  measure 
of  the  degree  of  intermetallic  mixing.  In  this 
graph  we  plot  the  average  number  of  nearest  neighbor  atoms  of  the  opposite  metal  for  atoms 
with  at  least  one  neighboring  atom  of  the  opposite  metal.  Since  the  timescales  of  our 
simulations  are  short,  only  a  very  small  amount  of  mixing  is  observed.  Cookoff  trajectories 
computed  for  nearly  a  nanosecond  did  not  show  a  significant  increase  in  reactivity  (see  Figure 
39).  Figure  46  shows  the  Ni  and  Al  surfaces  at  each  interface  for  the  final  shocked  state  of  the 
perpendicular  system.  The  layers  are  separated  in  the  figure  in  order  to  show  the  degree  of 
mixing.  What  is  most  interesting  about  these  results  is  that  the  planarity  of  the  second  interface 
is  disrupted  more  than  the  first.  We  believe  this  is  due  to  the  fact  the  shock  wave  moves  from 
the  material  with  higher  compressibility  into  the  material  with  lower  compressibility  at  interface 
2,  opposite  the  case  for  layer  1. 
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Fig.  44.  Shear  stress  and  temperature  profiles  as 
a  function  of  the  longitudinal  and  one  of  the 
transverse  dimensions.  The  snapshot  was  taken 
when  the  shock  wave  reached  the  free  surface. 
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Fig.  45.  Degree  of  mixing  of  each  metal  (Ni  and  Al)  with  the  other  metal.  The  quantity  plotted 
versus  time  is  the  average  number  of  nearest  neighbor  atoms  that  are  of  the  opposite  metal 
type.  Only  atoms  with  at  least  1  neighbor  of  the  opposite  metal  are  included  in  the  average. 


pertain  to  the  fully  shocked  state  of  the  system.  The  shock  traveled  from  the  left  to  right. 
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5.2  Toward  MD  Simulation  of  Shock-Induced  Chemical  Reactions  of  Titanium  and  Boron 


Molecular  dynamics  (MD)  simulations  using  a  reactive  force  field  were  used  to  study  the 
shock-induced  interfacial  chemistry  in  Ti/B  nanocomposites.  MD  trajectories  were  computed 
using  the  LAMMPS  software,  and  interatomic  forces  calculated  with  the  ReaxFF  reactive  force 
field.  Due  to  the  complexity  of  the  ReaxFF  potential  form  and  the  high  computational  cost  of 
computing  forces  with  this  force  field,  we  constrained  the  Ti/B  system  to  evolve  according  to 
the  Hugoniot  conservation  relations  using  either  the  NPzzHug  or  MSST  techniques.  The  usual 
method  for  computing  shock  properties  of  an  atomistic  model  system  is  to  perform  NEMD 
simulations.  The  drawback  of  this  technique  is  that  the  time  duration  of  the  simulation  is 
proportional  to  the  system  size.  As  a  result,  chemical  processes  that  occur  on  a  long  timescale 
cannot  be  modeled  with  standard  NEMD  techniques.  The  NPzzHug  and  MSST  methods  aim  to 
decouple  the  time  duration  of  the  simulation  from  the  size  of  the  system.  Both  methods  allow 
studies  of  the  chemical  reactivity  of  a  material  as  a  function  of  the  strength  of  mechanical 
impact  (external  pressure  in  the  case  of  NPzzHug,  and  shock  velocity  in  the  case  of  MSST).  The 
first  task  in  this  project  was  the  development  of  a  force  field. 

The  ReaxFF  force  field  parameters  were  fitted  to  atomic  coordinates,  atomic  charges, 
crystal  cell  parameters,  and  relative  potential  energies  taken  from  either  experimental 
measurements  or  high-fidelity  calculations  of  the  potential  energy  surface.  We  conducted  an 
extensive  set  of  quantum  mechanical  calculations  on  the  materials  involved  in  the  reaction  of 
titanium  and  boron.  These  calculations  include  single  point  energies  as  well  as  structural 
optimizations,  with  the  goal  of  including  the  following  in  the  training  set: 

1.  Fully  optimized  unit  cells 

2.  Unit  cell  compression  and  expansion 

3.  TiB2  crystal  vacancies 

4.  Gas  phase  PES  studies 

5.  Snapshots  of  a  periodic,  layered  Ti/B  system  undergoing  reaction 

Training  Set  Construction:  Solid  state  density  functional  theory  (DFT)  calculations  were 
conducted  with  the  ABINIT  plane  wave  code,  with  the  exception  of  those  pertaining  to  ooTi. 
This  structure  (for  unknown  reasons)  proved  difficult  to  model  with  plane  waves,  and  was 
therefore  modeled  with  localized  (Gaussian-type)  basis  sets  using  the  Crystal09  software.  The 
PBE  exchange-correlation  functional  was  used  for  all  calculations,  and  Trouiller-Martins-type 
pseudopotentials  were  used  for  all  plane  wave  calculations.  All  reported  energies  were 
computed  in  PI  symmetry  using  a  12x12x12  Monkhorst-Pack  (M-P)  grid  and  a  kinetic  energy 
cutoff  of  either  30  au  (for  supercell  structures)  or  60  au  (for  single  unit  cells).  Atomic  charges 
were  obtained  from  a  Bader  (atoms-in-molecules)  post-processing  analysis  of  the  electron 
density.  Single  molecule  studies  were  conducted  with  the  Gaussian03  software,  using  the 
LANL2DZ  basis  set  and  effective  core  potentials.  Periodic  Car-Parrinello  MD  simulations  were 
conducted  with  the  NWChem  pseudopotential  plane  wave  code  (£cutoff  =  20  au,  2x2x2  M-P  grid 
for  k  point  sampling). 
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Prior  to  the  work  reported  here,  there  were  no  ReaxFF  force  field  parameters  or 
parameter  training  set  for  the  reaction  of  Ti  and  B.  Parameters  existed  for  Ti  and  B  individually, 
but  not  for  their  interaction.  Much  of  the  current  effort  was  devoted  to  construction  of  a 
training  set  consisting  of  a  representative  set  of  atomic  configurations  and  their  associated 
energies  and  charge  distributions.  The  training  set  includes  the  following  bulk  Ti,  B,  TiB,  Ti2B, 
and  TiB2  structures: 

1.  Ti:  a  (hep),  (3  (bcc),  w(hcp)  phases,  as  well  as  simple  cubic  and  fee  structures 

2.  B:  the  rhombohedral  a  and  |3  phases,  the  high  pressure  orthorhombic  y  phase,  as  well 
as  simple  cubic,  fee,  bcc,  and  hep  structures 

3.  TiB:  orthorhombic,  bcc  and  hep  structures 

4.  Ti2B:  tetragonal  structure 

5.  TiB2:  the  highly  stable  hep  structure,  as  well  as  rhombohedral  and  tetragonal  structures 

The  most  stable  crystal  phases  at  room  temperature  and  standard  pressure  are  shown 
in  Figure  47.  The  high  pressure  phases  of  Ti  and  B  are  the  w  and  y  phases,  respectively.  For  all 
structures,  the  fully  optimized  unit  cell  structure  (coordinates  and  cell  parameters)  and  its 
energy  was  included  in  the  training  set.  The  cell  parameters  and  relative  energy  per  atom 
(fatom)  are  listed  for  Ti  and  B  in  Tables  V  and  VI,  respectively.  Tables  VII-IX  provide  the  cell 
parameters  and  formation  energies  (£»  for  the  TiB,  Ti2B  and  TiB2  stuctures.  Many  of  the 
structures  listed  above 
are  not  known  to  be 
stable  forms  of  the 
material  and  some 
were  found  to  have 
very  high  energies. 

These  structures  are 
included  in  the  training 
set  in  order  to  ensure 
the  force  field  will 
properly  predict  their 
instabilities.  The 

formation  energies 
calculated  for  hep  TiB2 
(-3.00  eV)  and 
orthorhombic  TiB  (- 
1.60  eV)  are  in  good 
agreement  with  the 
experimental  values  of 
-2.90  eV  and  -1.66  eV, 
respectively.50 


Fig.  47.  Most  stable  crystalline  structures  (at  ambient  temperature  and  pressure 
)  of  the  five  materials  included  in  the  training  set.  At  top,  from  left  to  right:  a 
2x2x2  supercell  of  aTi,  the  |3B  unit  cell.  At  bottom,  the  tetragonal  Ti2B  unit  cell, 
a  2x1x2  supercell  of  orthorhombic  TiB,  and  a  2x2x2  supercell  of  hep  TiB2. 
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Table  V  Crystal  phases  of  titanium. 


Phase 

Symmetry 

^atoms  /  UC 

0 

b 

c  Density 

fatom  (eV) 

a 

P63/mmc  (hep) 

2 

3.053 

4.524  4.368 

0 

p 

Im3m  (bcc) 

2 

3.349 

4.233 

0.098 

to* 

P6/mmm  (hep) 

3 

4.554 

2.808  4.730 

-0.010 

Pm3m  (cubic) 

1 

3.042 

2.824 

0.832 

Fm3m  (fee) 

4 

4.192 

4.316 

-0.033 

*  high  pressure  phase 

Table  VI  Crystal  phases  of  boron. 

Phase 

Symmetry 

N  atoms  / 

0 

b 

c  Density 

fatom  (gV) 

UC 

aa 

R3m  (rhomb.) 

12 

5.011 

2.506 

-0.134 

Pb 

R3m  (rhomb.) 

111 

10.139 

2.428 

0 

y* 

Pnnm  (orthorhombic) 

28 

5.025 

5.600  6.860  2.603 

-0.106 

Pm3m  (cubic) 

1 

2.067 

2.034 

1.582 

Fm3m  (fee) 

4 

2.919 

2.886 

1.257 

lm3m  (bcc) 

2 

2.128 

3.727 

1.883 

P63fmmc  (hep) 

2 

2.308 

2.805  2.774 

1.390 

a  symmetrized  crystal,  a  =  (3  =  y  =  58.557° 
b  symmetrized  crystal,  a  =  P  =  y  =  66.013° 
*  high  pressure  phase 


Table  VII  Crystal  phases  of  titanium  monoboride  (TiB). 

Symmetry 

^atoms  /  UC 

0 

b  c 

Density 

Ef  (eV) 

Pnma 

8 

4.577 

6.160  3.107 

4.451 

-1.604 

(orthorhombic) 

lm3m  (bcc) 

2 

2.833 

4.288 

0.040 

P63/mmc  (hep) 

2 

2.697 

4.315 

3.586 

1.225 

Table  VII  Crystal  phases  of  titanium  semiboride  (Ti2B). 

Symmetry 

Watoms  /  tic  a  b  c  Density 

Ef  (eV) 

14/mcm  (tetragonal) 

12  4.806  5.743  4.465 

-1.299 
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Table  VIII  Crystal  phases  of  titanium  diboride  (TiB2). 


Symmetry 

Watoms  /  UC 

a 

b  c 

Density 

ff(eV) 

P6/mmm  (hep) 

3 

3.045 

3.291 

4.366 

-2.997 

R3m  (rhomb.) 

9 

3.045 

11.610 

3.715 

-1.243 

P42/mnm  (rutile) 

6 

4.138 

4.635 

2.908 

5.758 

Compression/expansion  (c/e)  curves  were  computed  for  three  Ti  phases  (a,  3  and  ©), 
two  B  phases  (a  and  y)  and  hep  TiB2.  Cubic  and  uniaxial  compression  and  expansion  was 
computed  from  50%  compression  to  20%  expansion  in  increments  of  0.05.  Results  are  plotted 
in  Figure  48. 
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Fig.  48.  Compression/expansion  curves  for  Ti,  B  and  TiB2. 


Vacancy  energies  were  computed  for  removal  of  a  single  B  or  Ti  atom  from  2x2x2  and 
3x3x3  supercells  of  hep  TiB2.  In  these  calculations  the  atomic  coordinates  of  the  defective 
crystal  were  relaxed  with  the  cell  parameters  frozen  at  the  values  of  the  perfect  crystal.  The 
vacancy  energy  was  computed  according  to  the  following  equation 

^Ti/Bvac  =  [^TiBs(^‘l)  +  ^Ti/bI  “^TiBj(^) 

where  ETiB(N-l)  is  the  energy  of  the  defective  supercell,  EnB-M  the  energy  of  the  perfect 
supercell,  and  ETi/B  the  energy  per  atom  of  Ti  or  B  in  the  bulk  material  (a  phase  in  both  cases). 
Results  are  summarized  in  Table  IX. 


Table  IX  Titanium  diboride  (hep)  vacancy  energies. 


Supercell 

Vacancy  fraction 

^Ti  vac  (®V) 

ftwac  (eV) 

2x2x2 

0.042 

4.555 

3.130 

3x3x3 

0.012 

4.506 

2.952 

Although  the  MD  studies  we  planned  for  the  force  field  all  involve  condensed  phase  systems  it 
is  often  useful  to  include  work  on  gas  phase  systems  in  the  training  set.  The  main  advantages 
with  gas  phase  systems  are  they  are  often  smaller  and  hence  computationally  cheaper  to 
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model,  and  they  can  be  used  to  make  the  bond  stretch,  angle  and  torsion  terms  of  the  ReaxFF 
potential  form  more  straightforward  to  parameterize.  An  extensive  examination  of  gas  phase 
clusters  containing  only  Ti  and  B  yielded  several  stable  potential  energy  minima.  Unfortunately 
however  the  ground  electronic  states  of  these  Ti-B  clusters  are  rarely  singlet  spin,  making  their 
connection  to  the  condensed  phase  system  somewhat  complicated.  In  an  attempt  to  avoid  this 
difficulty,  we  investigated  a  Ti-B  cluster  containing  capping  H  atoms  (TiB4H8).  This  compound 
has  very  weak  Ti-B  bonding  however,  and  as  a  result  the  potential  energy  as  a  function  of  the  B- 
Ti-B  bond  angle  is  effectively  constant.  Ultimately,  we  decided  to  study  the  potential  energy 
surface  of  triplet  (ground  state)  Ti2B6,  and  include  these  results  (summarized  in  Figures  49  and 
50)  in  the  training  set.  This  compound  has  similarities  in  structure  to  the  TiB2  crystal:  six  boron 
atoms  in  a  hexagonal  ring  with  Ti  atoms  on  either  side  of  the  plane,  forming  a  line  that 
intersects  the  centroid  of  the  hexagon.  The  spin  density  resides  on  the  Ti  atoms.  The  B-B  and 
Ti-Ti  bonds  are  shorter  (by  0.16  and  0.19  A,  respectively)  and  stronger  than  those  of  the 
optimized  crystal,  because  the  atomic  coordination  numbers  are  significantly  lower  in  the 
cluster.  On  the  other  hand,  the  Ti-B  distances  (2.19  A)  are  very  similar  to  those  of  the  crystal 

(2.20  A). 


Fig.  49.  Potential  energy  surface  scans  of  the  triplet  Ti2B6  cluster  The  graph  at  the  top  shows  the  energy  as  a 
function  of  Ti-Ti  distance.  The  bottom  plot  graphs  the  energy  as  function  distance  of  the  boron  atom(s)  from  the 
midpoint  of  the  Ti-Ti  bond  vector,  which  is  the  center  of  mass  when  the  cluster  is  in  perfect  symmetry. 
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Fig.  50.  Potential  energy  surface  scans  of  the  triplet  Ti2B6  cluster.  The  energy  is  plotted  as  a  function  of  the  angle 
between  the  Ti  atom(s),  the  centroid  of  the  B  hexagon,  and  the  plane  of  the  B  atoms.  For  curve  a,  only  one  of  the 
Ti  atoms  is  displaced,  and  for  curves  b  and  c,  both  are  displaced  by  the  same  angle. 


Condensed.  ghase.  MD  sampling,:  Finally,  we  performed  a  high-temperature  Car-Parrinello  MD 
simulation  in  order  to  produce  amorphous  structures  to  use  for  both  force  field  training  and 
validation.  First,  a  layered  Ti/B  system  was  constructed  from  TA  unit  cells  of  aTi  (12  atoms)  and 
a  single  unit  cell  yB  (28  atoms).  The  cell  parameters  and  atomic  coordinates  of  the  system  were 
then  optimized  using  ABINIT  (2x2x2  M-P  grid,  20  au  energy  cutoff).  The  resulting  structure  is 
shown  in  the  far  left  of  Figure  51,  doubled  along  a,  b  and  c,  to  give  a  better  view.  During  the 
MD  simulation,  the  system  evolved  from  zero  temperature  to  3000  K  using  a  Nose-Hoover 
thermostat  with  a  time  constant  of  0.24  ps.  The  volume  of  the  simulation  cell  was  kept  fixed 
throughout.  Snapshots  of  the  system  along  the  trajectory  are  shown  in  Figure  51.  The  atomic 
coordinates,  cell  parameters,  relative  energies,  and  atomic  charges  for  each  of  these  snapshots 
were  included  in  the  training  set  for  force  field  parameterization. 


Fig.  51.  Snapshots  of  periodic,  layered  Ti/B  system  undergoing  thermally-induced  chemical  reaction.  The  original 
simulation  cell  (outlined  in  red,  green,  blue  and  black  lines)  was  doubled  in  each  dimension  for  visualization  clarity. 
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ReaxFF  Parameter  Optimization:  Training  set  consists  of  more  than  300  atomic  configurations, 
representing  both  single  molecule  and  condensed  phase  structures,  with  associated  energies 
and  charge  distributions.  Training  the  ReaxFF  potential  to  accurately  describe  the  solid  state 
chemistry  of  Ti  and  B  involves  optimization  of  78  force  field  parameters,  describing  Coulombic 
interactions  as  well  as  bond  stretching,  bending  and  torsion.  Parameter  optimization  is 
achieved  through  minimization  of  a  cost  function  using  a  single  parameter  search  algorithm.  In 
general,  an  adequately  parameterized  force  field  is  obtained  through  proper  description  of  the 
reference  properties  (the  training  set)  as  well  as  proper  weighting  of  these  properties  in  the 
cost  function.  For  example,  though  the  training  set  consists  of  approximately  300  atomic 
configurations  and  their  associated  properties  (e.g.,  potential  energies),  some  of  these 
configurations  (particularly  those  at  or  near  an  energy  minimum)  are  more  important  than 
others.  In  the  energy  cost  function,  the  structures  should  be  weighted  according  to  their 
"relative  importance."  Determining  the  relative  importance  of  each  data  in  the  training  set  is 
challenging,  and  unfortunately,  thus  far  our  optimization  attempts  have  yielded  unsatisfactory 
results.  Work  is  ongoing,  and  the  training  set  data  is  available  upon  request. 
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